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ABSTRACT 


The  problem  of  tracking  submarine  targets  with  passive 
sonobuoys  is  mathematically  modelled.  A  comprehensive  study 
is  made  of  all  the  information  available  in  the  acoustic 
signals  picked  uo  by  the  sonobuoys  and  of  the  usefulness  of 
this  information  in  the  estimation  orocess.  The  presence  of 
non 1 i near i t i es  in  the  tracking  model  leads  to  the  applica- 
tion of  nonlinear  estimation  theory,  Bayes  formulation  con- 
cepts are  applied  to  generate  approximate  solutions  and 
filtering  algorithms^  and  the  well  known  extended  Kalman 
filter  equations  and  higher  order  filtering  algorithms  are 
obtained  from  this  approach. 

The  concept  of  partitioning  the  measurements  is 
presented  and  shown  to  bring  advantages  in  computing  effi- 
ciency and  also#  for  nonlinear  measurement s f  in  tracking 
accuracy.  A  graphical  i nt eroret a t i on  of  the  action  of  Kalman 
filters  is  developed  and  provides  insight  into  the  impor- 
tance of  each  variable  in  the  filtering  process. 

Extensive  simulationsr  designed  to  test  the  performance 
of  the  algorithms  developed*  are  presented  in  graphical  form 
and  anal yzed. 
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SYMBOLS  AND  CONVENTIONS  FOR  TEXT 


1  -   Underlined  small  letters  represent  vectorsr  thusr 
x_   denotes  the  vector  x. 

2  -   Capital  letters  are    used  to  represent  matrices. 

3  -   A  caret  indicates  an   estimated   valuer   e.g./   x 
denotes  an  estimate  of  x.  • 


U       -   Time   points  are       represented   by   t,    or^   when 

k 

between   parentheses*   simply   by   k*  e.g.*  x_(  k )  denotes  the 

va  1  ue  of  x  at  t  i  me  t,  . 

—  k 

5   -   Estimates  specify  the  time  of  estimation  and   the 
amount   of   information  used*  e.g.#  x(k!j)  denotes  the  esti- 
mate of  X  at  time  t,  taking  into  account  all  the  measurements 
—  k 

UP  to  ana  including  t.  «  When  k  =  j  only  one  time  point  is 
indicated/  thus*  x(k)  aenotes  the  estimate  of  x  at  time  t, 
taking  into  account  all  measurements  up  to  and  including  t,  , 


6   -   Probability  density  functions  are    represented   by 

the   small   letter  p,     i.e.  p  i  (xly)  denotes  the  conditional 

x|y 

density  of  x  given  y.  when  no  confusion  is  possible/  this  is 
simplified  to  p(x|y). 


7   -   The  expectation  operator  is  represented  by  E  [   3; 
the  variance  by  Varf  1. 


a(k)     expected  value  of  the  estimated  state  vector  at 


\ 


a  (k)    random  pate  of  chanqe  of  target  soeed. 


a,(k)    random  rate  of  chanqe  of  target  heading. 
D         ■ 


a^(k)    random  rate  of  change  of  frequency  emitted 


b(k)     heading  of  target  at  time  t,  . 


b.(k)    bearing  measurement  by  buoy  i. 


b.(k)     second  moment  of  the  state  estimate  around   expected 
va 1 ue  a ( k  )  , 

c        average  speed  of  sound. 


d.(k)    distance  of  target  to  buoy  i. 

1 


kj 


Kronecker  delta;  =lifk=j;=Oifk=j 


e(k)     estimation  error  vector. 


e(k)     expected  value  of  the  estimation  error 


f  (k)    rest  frequency  emitted  by  the  target 
o 


f. (k)    frequency  measured  by  buoy  i 

X 


G ( k  3     gai  n  mat  r i  x  , 


H(k)     observation  matrix. 


compression/expansion  factor 


P(kJ     estimation  error  covariance  matrix. 


Q(k)     state  excitation  covariance  matrix. 


R(k)     measurement  noise  cavariance  matrix. 


s  (  k  )     speed  of  target. 


s,(k)    relative  velocity  of  target  towards  buoy  i 

1 


a  (k)    standard  deviation  of  a   (k), 
s  s 


a,  (k)    standard  deviation  of  a^  (k). 

D  D 


aAk)         standard  deviation  of  a   (k). 


t  i  me  del  ay  . 


T.       time  delay  difference  between  the   signals   received 


by  buoys  i  and  j . 


time  spent  in  prediction. 


time  soent  in  estimation. 


v(k)     vector  of  q  random  measurement  noise  signals 


wCkJ     vector  of  m  random  forcing  inputs 


xCk)     state  vector  of  dimension  n. 


x(k)     x-position  of  target. 


x.(k)    x-Dosition  of  buoy  i. 

1 


x(k)     estimated  state  vectoi 


y(kj      [n (n+3 ) /2] -d i mens i ona 1  vector  containing   the   state 
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variables  and  the  distinct  conponents  of  the  matrix  P. 


y(k)     y-position  of  target. 


y. (k)    y-position  of  buoy  i. 


z(k)'    vector  of  q  measurements. 


Z        set  of  all  measurements  up  to  and  including  time  t  t 

k 

i.e.   tz(U),  z(l)r  ...f  zCk)l  . 
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I.   INTRODUCTION 


This  work  was  oriented  to  the  problem  of  optimally 
estimating  characteristics  and/or  parameters  of  a  certain 
class  of  nonlinear*  dynamic/  stochastic  systems  from 
observed  output  sequences  which  are  noise  corrupted  non- 
linear functions  of  the  system  states. 

In  particular*  attention  was  directed  to  the  problem  of 
real-time*  short-range  optimal  localization  and  tracking  of 
a  submarine  target  by  passive  acoustical  means.  The  most 
accurate  and  reliable  estimates  of  the  target's  parameters 
(position*  heading*  sceed)  are  sought  by  optimally  process- 
ing the  measurements  obtained  from  the  acoustic  signals 
transmitted  by  the  target  itself  and  received  by  special 
sonobuoys . 

In  the  unclassified  literature  there  are  presently 
available  methods  of  passive  tracking  of  submarine  targets. 
Reference  [IJ  utilizes  doool er-sh i f t ed  frequency  measure- 
ments obtained  from  a  group  of  sensors*  12]  uses  bearing  and 
dopp 1 er-sh i f t ed  frequency  measurements  obtained  from  one  or 
more  sonobuoys.  These  methods  provide  very  good  first 
approximations  to  the  solution  of  the  problem  and  many  of 
their  concepts  and  notation  are    used  in  this  work. 
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As  a  first  step*  an  attemot  was  made  to  produce  a 
comprehensive  study  of  all  the  information  available  in  the 
acoustic  signals  picked  up  by  the  sonobuoys  and  its  useful- 
ness  in  the  estimation  process. 

The  estimation  algorithms  were  developed  taking  into 
account  the  advantages  of  processing  information  as  soon  as 
it  becomes  available.  Flexibility  to  move^  remove  and 
include  sensors  and  measurements  at  any  time  was  obtained. 
Constraints  were  included  to  account  for  limitations  of 
practical^  inexpensive  and  presently  available  sonobuoys  and 
computing  eguipment. 

In  Chapter  II  the  problem  is  described  in  mathematical 
terms  using  state  space  techniques  and  the  initial  assump- 
tions and  constraints  are  given.  The  model  used  for  the 
target  is  similar  to  the  model  presented  in  [21 r  but  with  a 
different  set  of  states.  The  measurements  obtainea  from  the 
signals  received  by  the  sonobuoys  are  characterized  and 
their  functional  relationships  to  the  parameters  of  the  tar- 
get are    anal y zed . 

Chapter  III  discusses  the  general  theoretical  solution 
for  the  problem  and  the  difficulties  of  implementing  it. 
Approximate  solutions  are  then  sought  and  processing  equa- 
tions are  developed.  The  special  vectors*  matrices  and 
relations  characteristic  of  the  tracking  problem  are 
prepared  for  the  application  of  the  filtering  eauations. 
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In  Chapter  IV  the  concept  of  partitioned  measurements 
is  introduced  and  analyzed.  The  advantages  in  accuracy/  com- 
puting efficiency  and  filter  flexibility  are  stressed. 

Practical  and  graphical  interpretations  of  the  estima- 
tion process  help  in  visualizing  nonlinear  problems  and 
motivate  the  adoption  of  iterative  techniques.  Increased 
accuracy  and  possible  convergence  improvements  are  shown  to 
result  from  this  approach. 

Chapter  V  introduces  the  interactive  computer  program 
that  allowed  the  application  of  the  ideas  described  in  the 
previous  chapters  to  the  tracking  problem.  Selected  simula- 
tions are    discussed  and  presented  in  graphical  form. 

This  research  covers  the  subject  in  a  much  more 
comprehensive  way  than  earlier  works/  [11/  12],  [31.  The 
results  obtained  by  using  the  methods  and  ideas  collected 
and  developed  in  this  study  show  many  ways  of  improving  the 
tracking  accuracy/  of  increasing  computing  efficiency/  of 
reducing  divergence  problems  and  of  obtaining  faster  conver- 
gence. 

The  importance  of  having  the  buovs  placed  in  adequate 
positions/  of  havinq  an  adequate  model  for  the  target/  and 
of  having  frequent/  varied  and  precise  measurements  is 
st  ressed . 


15 


The  final  chapter  summarizes  the  results  of  this  inves- 
tigation and  presents  the  conclusions  and  suggestions  for 
f urt  her  study . 
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II.   PROBLEM  FORMULATION 


A.   THE  ESTIMATION  PROBLEM 

A  general  nonlinear^  stochasticr  dynamic  system  observed 
by  a  set  of  nonlinearr  nonst at i onary  measurements  can  be 
represented^  for  our  purposes^  in  state-space  discrete  form 
by  t  he  equat  i  ons 


x(k  +  l)  =  f(x(k),w(k),t,  ,t,  ,    ) 


_z(k)  =  jiC)<_(k),v^(k)  ,t  ) 


(2.1) 


(2,2) 


where   x_(  k )  is  the  state  vector  of  dimension  n^ 

w^(k)  is  the  vector  of  the  n  random  forcing  inputs^ 

_z  ( k )  is  the  measurement  vector  of  dimension  Or 

v^(k)  is  the  vector  of  p  random  measurement  noises* 

_f(  )  and  h(  )  are  general  vector  functions. 

The  general  estimation  problem  consists  of 

---  knowing  the  statistical  characteristics  of  w(k)  and 

j^CkJ   and  having  a  statistical  representation  of  the  initial 

conditions  x_(  0 )  of  the  system* 

---  processing  the  observations  z^(.0)  ,     z_(  1  )  *   .../   £U) 

in   a  statistically  optimal  way  to  obtain  the  best  estimate* 

in  some  sense*  of  the  system  state  x(k)  at  time  t,  . 

—  k 

If  k  >  j  the  Drocessing  is  called  prediction*  if  k   <   j 
it  is  called  smoothing  and  if  k  =  j  it  is  called  filtering. 
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In  this  work  some  departures  from  the  completely  general 
nonlinear  case  were  taken.  The  first  assumption  is  that  the 
measurement  noise  is  additive  and  that  the  system  equations 
are  linear  with  respect  to  the  random  input  vector  w_.  Equa- 
tions id,\)    and  (2,2)  can  then  be  written  as 


x(k  +  l)  =  f  (x(k),t,  /t^,  )  +  g(x(k),t  ,t    ).w(k) 

—  k,   lC+1        ^     —  1,     1,4.1     — 


k   k+1 


£(k  +  l)  =  tl(jl(k),t  )  +  v_(k) 


(2.3) 


(2. a) 


where  f_(  J  and  h^(  )  represent  general  vector   functions   and 
g(  )  a  general  matrix  of  functions. 

Secondly^  the  measurement  noise  and  random  forcing  input 
are  assumed  to  be  uncorrelated  in  timer  zero  mean,  discrete 
Gaussian  sequences,  independent  of  one  another  and  indepen- 
dent of  the  initial  condition  of  the  state  of  the  system, 
i.e.. 


E  [w(kn  =  0 


E [v(k)]  =  0 


(2.5) 


E  lwlk)w^(  j  )]  =  Q(k).  6, 


E  [v(k)v^(  j  )]  =  R(k).  6, 


Elw(k)v^(j)]  =  0     E(x(0)w^(k)]  =  0     E[v(k)x^(0)]  =  0 


kj 


1  i  f  k  =  j 

0  i  f  k  ?i  j 
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B.   THE  PASSIVE  TRACKING  PROBLEM 

Many  practical  situations  can  be  formulated  as  estima- 
tion problems  and  characterized  by  the  equations  described 
above.  The  situation  that  motivated  this  study  assumes  a 
submarine  target  following  (most  of  the  time)  an  approxi- 
mately constant  speed  and  constant  heading  path  in  a  field 
of  passive  sonobuoys. 

The  target  unintentionally  emits  acoustic  signals  which 
are  picked  up  by  hydrophones.  The  transduced  signals  are 
sent  UP  by  wire  to  the  buoys  and  then  frequency-modulate  VHF 
carriers  which  are  transmitted  by  the  buoys  and  received  at 
a  nearby  ship  or  aircraft.  After  the  recovery  of  the  signals 
data  processing  equipment  generates  measurement  values  which 
contain  numerical  information  about  the  target  parameters. 

The  first  step  in  the  analysis  of  the  signals  collected 
by  the  sonobuoys  is  the  attempt  to  detect  the  presence  of  a 
real  target • 

After  a  target  is  detected^  frequency  spectrum  informa- 
tion is  added  to  intelligence  data  in  an  attempt  to  classify 
it.  Other  measurements  and  information  from  any  other  source 
are  also  used  to  obtain  a  first  approximation  to  the 
target's  parameters  (position^  headingr  speed)^  generally 
through  the  use  of  least  mean  square  techniques  [11 f  [3]. 

The  third  phase  of  the  processr  if  a  wartime  condition 
doesn't   existf   is  the  track i no  of  the  tarqet.  That  iS/  the 
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determination  of  its  path  while  reducing/  if  possible/  the 
initial  uncertainties  about  its  parameters/  in  real-time/ 
through  the  optimal  processing  of  the  sparse  information 
provided  by  the  sonobuoys.  The  passive  characteristic  causes 
severe  limitations  but  is  necessary  to  avoid  revealing  to 
the  target  the  fact  that  it  is  being  observed. 

This  investigation  was  devoted  to  the  tracking  phase  and 
it  was  assumed  that  the  main  characteristics  of  the  buoys  as 
well  as  their  positions  are  known, 

C.   THE  MODEL  OF  THE  TARGET 

The  problem  is/  o^  course/  three-dimensional.  Neverthe- 
less the  measurements  are  not  a  direct  function  of  the 
target's  deoth  but  of  the  difference  in  depth  between  the 
target  and  the  hydrophones.  Anticipating  the  accuracy  of  the 
measurements  and  the  precision  of  the  data  processing  eauip- 
ment  and  algorithms/  it  is  justifiable  to  consider  only  two 
dimensions  initially/  adding  depth  as  the  third  dimension 
whenever  necessary  ana  with  no  conceptual  aifficulty. 

The  target's  path  is  freguently  one  with  nearly  constant 
soeed  and  heading/  disturbed  by  currents  and  random  thrust 
and  control  variations.  Intentional  maneuvers  whicn  have  no 
evasive  purposes  are    normally  simple  and  smooth. 

A  basic  plant  havino  as  states  two  position  coordinates 
(x  and  y)/  heading  (b)  and  soeed  (s)  of  the  target/  describ- 
ing a  constant  velocity  path/   as   shown   in   Fig   1/   seems 
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adequat  e . 

Random  forcing  inputs  should  be  considered  to  account 
for  noise  processes  and  imperfections  of  the  models  such  as 
target  maneuvers. 

As  suggested  by  12J  these  random  forcing   functions   can 

be   approximately   represented  by  two  independent  zero-mean^ 

pi  ecew  i  se-const  ant  random  rates  of  change/  a   si^d  a  >       act- 

s        b 

i ng  on  the  speed  and  heading  of  the  target. 

The  formulation  of  this  basic  plant  in  discrete  state- 
space  formf  similar  to  Equation  (2.3)/  can  be  obtained  in 
the  following  way 


'k+1 


sCk+l )  -  s(k)  = 
b(k+l)  -  b(k)  = 
x(k  +  l  )  -  x(k)  = 


a  (t)  dt  =  (t,  ,   -  t  )  .  a  (k) 

s  k+1     k    s 


(t 


k+1 


a,  (t)  dt  =  (t.     -  t  ).a  (k) 
b  k+1      k    b 


(t 


k+1 


s  (t )  cos  b( t )  dt  = 


'k+l 


[s(k)  f  (t  -  t,  )a  (k)]  cos  [b(k)  + 
k   s 


+  (t  -  t,  )a,  (k)]  dt  = 
k   b 


"k+1 


s(k)  cos  b(k)  dt  + 


ft 


k+1 


^k+1 


(t  -  t,  )a  (k)  cos  b(k )  dt  - 
k   s 


(t  -  t,  )a,(k)  s(k)  sin  b(k)  dt  + 
k   b 


tk 


=  (t,  , ,   -  t,  )  s  ( k)  cos  b  ( k)  + 

k+1      k 
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target  position 
at  time  tj^ 


\ 


A 

t  buoy    j 


\ 


"tjuoy    i 


\ 


I  \ 

I 


x.(k; 

1 


x.(k) 
J 


X 


Figure   1    :    Geometry   of   target  and  sensors 
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+  ^   (t 

2    k+1 


-  t  )^  [  a  (k)  cos  b(k) 

k      s 


-  a  (k)  s(k)  sin  b(k)J  +  HOT 

b 


(2.6) 


where  HOT  represents  Higher-Order  Terms. 


In  the  same  way^ 


rt 


y(kfl)  -  y(k)  = 


k+1 


s(t)  sin  b(t)  dt  = 


=  (t,  ,  -  t,  )  s(k)  sin  b(k)  + 

k+1     k 


+  -I  (t     -  t  )^  [a  (k)  sin  b(k)  + 
2     k+1     k      s 


+  a,(k)  s(k)  cos  b(k)]  +  HOT 

D 


(.^.7) 


If  Ct.^T  -  t.  )a  (k)  and  (t,  ,,   -   t,  )a,  (k)  are       suffi- 

k+1     k   8  k+1       k    b 

ciently   small   so   that  the  higher  order  terms  in  Equations 
(2.6)  and  (2.7)  can  be  neglected*  one  has*  in  vector  form 

x(k  +  l)  =  f(xCk),t  ,t,^,  )  ■»•  g(x(k),t.t,   ).w(k) 
—         —  —     k   k+1      k   k+1   — 


where 


x(k)  = 


x(k) 
y(k) 
s(k) 
b(k) 


(2.8) 


w(k)  = 


a  (k) 

s 


a   (k) 

b 


(2.9) 
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f(x(k),t    ,t     )  = 

-  -      k+1   k 


(2.10) 


x(k)  ♦  (t     -  t  )  s(k)  cos  b(k) 

k+1    k 

y(k)  +  (t     -  t  )  s(k)  sin  b(k) 

k+1    k 

s(k) 
b(k) 


and 


g(x(k),t    ,t  )  = 
—  —     k+1   k 


1  2 

"^  ( t  ,  , ,  -  t,  )  cos  b  (  k  ) 

2  k+1     k 

i^^k+1  -  \^'^^"  ^^^^ 
^^k+i  ■  ^k^ 


(2.11) 


-^(t,.,  -t,  )^s(k)  sinb(k) 
2   k+1     k 


^(t  ,  ,,  -  t,  )^  s(k)  cos  b(k) 

2   k+1     k 

0 


^^k+1  -  ^k^ 


When   frequency   measurements   are    to    be    processed 

directly^   as  shown  in  Section  IlfD^l^D^  an  additional  state 

variable  must  be  included  ---  the  rest  frequency  f    emitted 

o 

by  t  he  target  . 


This  frequency  is  assumed  to  be   approximately   constant 

with   a   small   random   disturbance  a  ^r  that  is  assumed  to  be 

zero-mean^  p i ecew i se-const ant »  and  independent  of  a  and  a.  . 

s        D 

The  expanded  state  variable  formulation  is  given  by 


x'^Ck)  = 


x(k) 


f  (k) 
o 


(2.12) 


2a 


w^Ck)    = 


w(k) 


a^(k) 


(2.13) 


-    <e  r    e 


x^  (k  +  n    =    f*==  (x*==(k),t,_,,    ,t,   )     + 
-  —      —  k+1       k 


+    q^(x^  (k)rt  ,^wt,  ).w®(k) 
^      —  k+1      k      — 


(2.1^) 


where 


f  ^(x^(k),t,^T   't,  )    = 
—       —  k+1       k 


k+1        k 


f     (kJ 
o 


(2.15) 


/(x.^k),t^^^,t^)    = 


S^^^^^'Vl'^^ 


^^k+1-    \' 


(2.16) 

D.   MEASUREMENTS 

The  songbuoys  can  perform  one  or  both  of  the  following 
tasks : 

---  pick  up  underwater  sound  signals 

---  indicate  the  approximate  direction  of  the  vectorial 
sum  of  the  received  acoustic  signals. 
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The  target  emits  a  signal  r(t)  which   is   distorted  and 

modulated   by   the  propagation  medium  between  the  source  and 

the  hydrophonef  and  by  extraneous  sound  sources   present  in 
the  ocean. 

When  all  of  these  disturbances  are  of  relatively  mild 
strengthr  the  signal  received  by  a  stationary  buoy  i  will  be 
approximately  an  attenuated  and  delayed  version  of  the  pre- 
viously emitted  signals  i.e./ 


r  (t )  =  a  (t).r(t  -  t  (t)) 
i        i  i 


(2.17) 


where  a  (t)  is  an  attenuation  factor  and  x. (t)  is   the   time 
i  1 

de 1  ay  gi  ven  by 


T  (t)  =  d  (t  -  T  (t))  /  c 
i         i       i 


(2.18) 


that  \Sf     the  distance  d   between  the  target  and  the  buoy   at 

i 

the   time   of   emission   divided   Dy  the  average  velocity  of 
sound  in  the  medium. 

Since  the  target  velocity  is  much  smaller  than  the  sound 
velocity/  for  small  distances  and  delays  one  has 


d.(t  -  T  (t))  =  d  (t)  +  s  (t).x  (t)  (2.19) 

1       i         i        i      i 

where  s  (t)  is  the  relative  speed  of  the  target  toward   buoy 
i 


i  as  shown  in  Fig  2 . 


ThuS/  from  Equations  (2.18)  and  (2.19)  one  has 


T.  (t )  =  d. (t )  /  (c  -  s  (t)  ) 
IX  i 


(2.20) 
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position  of  target 
at  time  t-X,(t)  . 


position  of  target 
at  time  t . 


buoy  j 


uoy  1 


\ 


d(t-t(t)) 
i    i 


■^ 


Figure  2  :  Target-sensors  geometry  -  distances 
and  relative  velocities. 
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and  the  received  signal  is  qiven  bv 


p.(t)  =  a.(t).r[t  -  d.(t)/(c  -  s  (t))]       (2.21) 

IX  1  1 

If  this  signal  is  recorded  from  t  =  t   to  t  =   t    +   T^ 

o  o 

one  has  recorded 


r(t   +t')=a(t   tt').rlt   +t'- 
10  i   o  o 


-d  (t   +  t  •  )/(c  -  s  (t   +  t •))] 
i   o  i   o 

0  <  t  •  <  T 

If  during   this   oeriod   T   the   relative   target   speed 
towards  buoy  i  doesn't  change  significantly^ 


s(t       +t')=s(t) 
i      o  i      o 


d.  (t       -^    t  ')    =    d    (t     )    -    t  •  .s     (t     ) 
1      o  i       o  i       o 


and    hence 

r.(t       tt')=a.(t  +t').rrt       -d(t)/(c-s(t)J 
1       o                               i       o  L.O  i       o  i      o 

+    t'.[l     +  s.(t     )/(c    -    s     (t     ))]] 


1        o 


1        o 


or 


r.(t       +    t*)    =    a.(t       +    t').rtt       -    d    (t    )/(c    -    s    (t     ))    + 
io  io  oio  io 


+    c.t  '/(c    -    s     (t     ))] 

i        o 


(2.2?) 


In  tquation  (.2,22),  the  second  term  inside  the  brackets 
is  a  fixed  time  delay.  The  third  term  is  a  variable  time 
delay  and  can  be  seen  as  a  comoressi on/expans i on  term,  show- 
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ing   the  different  value  of  the  variable  time  for  the  origi- 
nal and  the  received  signals. 

If  the  emitted  signal  had  a  single  frequency  component  w 
one  wou Id  have 


r(t       ■♦•t')=a(t       ftM.cosCwCt       +t')) 
i-     o  i       o  i       o 


where 


cos(w    (t       +    t*))    =    cos(4)        +    w    .t')    = 
i       o  io  i 


cosfw    .t       -    w    .d    (t     )/(c    -    s    (t    ))    + 

CO  o       i      o  i       o 

•♦•    w    .c.  t  •  /(c    -    s     (t     )  )]     = 

o  i       o 


C0S((f)        +    w    .c.t'/(c    -    s     (t     ))] 
io  o  i       o 


and  the  received  frequency  would  be 


w   =w.c/(c-s) 
i     o  i 


f   =  f  .c  /  (c  -  s  ) 
i     o  i 


(2.23) 


The  signals  collected  by  various  sonobuoys  are  sent  to  a 
ship  or  aircraft  normally  equipped  with  equipment  and  skills 
to  extract  the  information  described  below. 

I  .   I  so  1  at ed  Buoys 


By  processing  the  signals   collected   by   each   buoy 
alone*  one  can  obtain  the  following  noisy  measurements; 
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a.   Bearing  Indication 

Bearing  is  given  directly  by  the  buoys  or 
preprocessed  (filtered)  to  reduce  noise  effects. 

Since  the  bearing  is  obtained  from  a  vectorial 
sunif  it  can  be  reasonably  accurate  only  when  the  acoustic 
noise  is  approximately  omnidirectional  or  much  weaker  than 
the  target  signal.  For  some  sensors  this  vectorial  sum  can 
be  limited  to  a  selectable  frequency  band#  thereby  easing 
the  noise  problem. 

The  relationships  between  a  bearing  indication 
obtained  by  a  buoy  i  and  the  states  of  the  target  is  shown 
in  Fig  i    and  given  by 

b.(k)  =  arctanflyCk)  -  y  (k)]/[x(k)  -  x  (k)]"]  +  v(k) 
X  *-  i  i    -* 

(?.2a) 

where  v(k)  must  account  for  all  the  noise  aisturbinq  this 
measurements  including  acoustic  noise/  transducer  inaccura- 
cies* preprocessing  noise  and  errors  due  to  the  finite  speed 
of  the  sound  in  the  water.  This  last  factor  is  caused  by  the 
time  delay  between  the  emission  of  a  sound  by  the  target  and 
reception  at  the  hydrophone.  At  the  end  of  this  period  the 
target  and  the  hydrophone  have  slightly  different  relative 
pos  i  t  i  ons , 

According  to  [2]  typical  accuracies  of  +  5 
degrees  are  common  for  strong  sianals  and  inexpensive  DIFAR 
buoyS/  and  as  cited  by  13]  inaccuracies  down  to  +   1   degree 
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y(k) 


y.(k) 

1 


b^(k)  -  b(k) 


s.(k)  ^ 


b(k) 


x(k)  x(k) 
i 


X 


Figure  3  :  Bearing  measurement, 
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OP   less   can   be   obtained   with   arrays  of  hydrophones  and 
proper  preprocessing. 

b.   Frequency  Indication 

By  analyzing  records  of  length  T  seconds  of  the 
acoustic  signals  picked  up  by  the  sonobuoys  with  Fast 
Fourier  Transform  algorithms/  one  can  obtain  approximate 
frequency  spectrum  descriptions  with  a  resolution  of  1/T 
Hertz. 

By  detecting*  measuring  and  possibly  tracking 
[4]  some  of  the  strongest  frequencies  contained  in  the  sig- 
nal* very  useful  data  about  the  target  parameters  is 
obtained  since  the  received  frequencies  are  functions  of  the 
originally  emitted  frequencies*  the  speed*  the  heading  and 
the  position  of  the  target*  as  shown  in  the  relationship 


f  (k)  =  f  (k).c  /  (c  -  s  (k))   +   v(k) 
1        o  i 


(2.^5) 


where 


s  .(k)  =  -s(k) .cos [b(k)  -  b  (k)l 
1  i 


(2.26) 


is  the  relative  velocity  toward  buoy  i*  as-'Shown  in  Fig  3. 

The  term  v(k)  now  has  to  account  also  for  the 
errors  introduced  in  the  various  paths  of  the  signal  from 
the  target  to  the  ship  or  aircraft*  and  for  variations  in 
t arget -sensor  geometries  during  the  record  time  T. 
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As  suqgested  in  []],[?]       and   (33 r   inaccuracies 

ranging   between   0,5  and  0.01  Hertz  can  be  obtained  in  most 

situations  for  values  of  f   of   the   order   of   hundreds   of 

o 

Hertz. 

c.   Signal  Strength  Indication 

The  strength  of  the  signals  varies  not  only  with 
the  strength  of  the  originally  emitted  sound  but  also  with 
the  distance  between  the  target  and  the  buoy.  The  signal 
strength  also  depends  on  the  random  aspects  and  the  direc- 
tivities of  the  targets  the  hydrophonesr  and  of  the 
transmitting  and  receiving  VHP  antennas.  Random  absorotion 
and  scattering  in  the  various  media  from  the  target  to  the 
processing  equipments  may  also  have  a  significant  influence 
on  the  signal  presented  to  the  processors. 

A  mathematical  formulation  of  all  these  effects 
would  be  extremely  difficult  task  and  would  obscure  the 
desired  distance  information.  For  this  reason  the  informa- 
tion contained  in  the  strength  of  the  collected  signalsr 
which  can  be  used  by  an  experienced  and  ski  Ilea  audio  opera- 
tor,  was  not  incorporated  in  this  study. 

Attenuation  factors  are  thrown  away  by  setting 
the  strength  of  the  received  signals  at  a  good  working 
1  eve  I  . 
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2.   Two  or  More  Buoys 

Considering  the  signals  received  by  two  or  more 
buoys^  indications  can  be  obtained  which  may  be  used  in 
addition  tOr  or  instead  of/  the  measurements  available  from 
a  s  i  ngl e  buoy . 

a.   Frequency  Difference  Indication 

To  process  frequency  measurements  obtained  from 
one  buoy  requires  inclusion  of  the  rest  frequency  as  one  of 
the  states.  This  is  the  case  because  this  measurement  is 
very  sensitive  to  errors  in  the  assumed  rest  frequency^ 
i  .e.  f 


f  .  =  f  .c  /  (c  -  s.  ) 

1      o  1 


Af 


Af.c/(C-S.)=   Af 
O  1  o 


Having  the  frequencies  received  by  two  buoys  and 
taking  their  difference  one  gets 


f.  -  f.  =  f  .c/(c  -  s.)  -  f  .c/(c  -  s.)  = 

1      J     o  1       o  j 


f  .c.(s.  -  s.)  /  rOc  -  s.).(c  -  s.))  = 

013         1       J 


f  . (s  .  -  s. )/c 

o    1    J 


(2.27) 


The  sensitivity  of  this  difference  to  errors   in 
the  assumption  of  t^e    rest  frequency  is  given  by 


A(f.  -  f.  )  =   Af  .  (s  .  -  s  )/c 
1     J         o    1     j 


3a 


or  is  reduced  by  a  very  large  factor  since  c  is  of  the  order 
of  1500  m/s  and  (s.  -  s.)  is  normally  very  small. 

This  way/  if  one  has  a  reasonable   approximation 

for   f  t       instead  of  processing  two  measurements  in  a  five- 
o 

state  plant  one  can  process  only  one  f reguenc y-di f f erence 
measurement  in  a  four-state  olant  with  greatly  reduced  com- 
puting time  and  hopefully  not  a  detectable  reduction  in 
accuracy • 

The  variance  of  the  errors  in  each  measurement 
must  be  added  to  give  the  variance  of  the  error  in  the  fre- 
quency difference  measurement. 

b.   Frequency  Ratio  Indication 

If  one  divides  the  frequencies  measured  at  two 
sonobuoysr  a  new  relationship  is  obtained  that  retains  the 
information  on  position/  heading  and  speed  of  the  target/ 
but. is  independent  of  the  rest  freguency  emitted 


f.  /  f.=  [f  .c/(c  -  s.)l.[(c  -  s.)/(f  .c)l  = 
3-      J      o  1  JO 


=  (c  -  s  )  /  (c  -  s  ) 
J  i 


(2.28) 


or 


f./  f.=  I  +  (s.  -  s.)/c  .  [1  +  s  /c  +  (s  /c)^  +  ...1 

1     J  1     J  1        i 


=  1  +  (s   -  s  )  /  c 
1     J 


(2.29) 
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As  with  the  frequency  difference  indication^ 
this  measurement  can  be  used  in  a  four-state  plant  instead 
of  using  two  frequency  measurements  in  a  five-state  plant. 

If  this  ratio  is  generated  by  the  division  of 
the  frequencies  obtained  from  the  analysis  of  each  of  the 
signals/  the  characteristics  of  the  frequency-ratio  measure- 
ment noise  is  complicated  and  state  dependent.  If  however 
this  relation  is  obtained  as  described  in  the  following  dis- 
cussion, this  uncomfortable  situation  can  be  avoided. 

c.   Time  Delay  Indication 

As  shown  at  the  beginning  of  this  section  on 
measurement s f  the  signal  collected  by  a  sonobuoy  is  a  noisy, 
delayed,  attenuated  and  compressed/expanded  reproduction  of 
the  original  signal  emitted  by  the  target.  The  delay  is  due 
to  the  finite  time  the  signal  takes  to  reach  the  buoy  and 
the  compression/expansion  is  caused  by  the  variation  of  this 
delay  with  time  (doppler  effect). 

The  signals  received  by  two  sonobuoys  have  dif- 
ferent noise  contributions  and  different  delay  and 
compression/expansion  factors  even  after  their  strengths  are 
equal i  zed. 

Supoose  one  recorded   the   signals   received   oy 

buoys   i   and   j   from   t=ttot=t   +T,  The  signals  are 

o  o 

related  to  the  original  signal  as  shown  in  Equation   (£.^2), 
i.e., 
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p.(t   ♦  t')  =  a,(t   +  t').r[t   -  d.(t  )/(c  -  s  (t  ))  + 
lo  io  oio  io 


t  c.f /(c  -  s.  (t  ))] 

X    o 


(2.30) 


r.(t   •»■  t')  =  a.(t   +  t').rtt   -  d  .  ( t  )/(c  -  s  (t  ))  + 
JO  JO  ojo  JO 


+  c.t  '/(c  -  s. (t  ))] 

J   o 


(2.31) 


If  one  now  amolifies  these  two  signals  to  a  good 
working  level/  correlate  them  by  evaluating 


(t)  = 


P  '.  T 


r.(t   +  t').r.(t   +  t'  +t) .dt ' 

1    o  JO 


(2.32) 


and  look  for  the  value  of   x  that  maximizes  p  (x)  one   finds 
that : 


(1)  if  s.  (t  )  =  s  (t  )  andf    consequently/  the  fre- 

1    o        JO 

quency   shifts  &re    aoout  the  same  in  both  signals^  then  r  (t 

i   o 

+  t')  will  be  simolv   an   aoproximate   delayed   or   advanced 

replica   of  r  (t   +  t')r     as  can  be  seen  from  the  above  equa- 
j   o 

t  i  ons . 

In  this  case  the  maximum  of   p  ( x)  is  easy   to 

find  and  corresponds  to  the  value  x      at  which 

max 

t      -    d.(t    )/(c    -    s.  (t    ))    +    c.(f    fx        )/(c    -    s    (t    ))    = 

o  1     o  i       o  max  i     o 

=    t       -    d    (t     )/(c    -    s     (t     ))     +    c.t'/(c    -    s     (t     )) 

ojo  jo  jo 


solving    for    x  yields 

max 


X         =     [d    (t     )     -    d    (t     )]     /    c 

max  10  jo 


(2.33) 
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(2)     if    s.(t     )     i    s.(t     ) ,       n(T    )       does       not       have       a 

1    O         JO       ^ 

simply-determined  maximum  value  because  the  two  signals  have 
frequency  components  whose  phases  are  shifted  by  different 
amount  s . 

If  ones  now  change  the  time  scale  of   one   of 
the  signals*  say  r.(t   +  t*)#  one  has 

JO 

r'(f)  =  r.(t   +  m.t')  = 

J  JO 

=  a.-rfct   -  d.(t  )/(c  -  s.(t  ))]  + 
J    L.  o     JO  JO 

+C.m,t'/(c-s.(t))l  (2.3^) 

J   o  -' 

Comparing  this  equation  with  Equation   (2,30) 
one  sees  that  if  m  can  be  found  such  that 


m.c/(c  -  s  .(t  ))  =  c/(c  -  s  (t  )) 
JO  i   o 


or 


m  =  [c  -  s  .(t  )]  /  [c  -  s  (t  )]  (2.35) 

JO  ID 


one  again  has  case  (1)  and   the   correlation   function   will 
give  a  maximum  at  approximately 


T    =  Cd .(t  )  -  d  .(t  )1  /  c  (2.36) 

max     1   o      JO 


It  is  interesting  to  note  that   Eauation   (2.35) 

is   the   same  as  Eauation  (2.28).  The  compression/expansion 

factor  is  thus  the  ratio  between  the  frequencies  received  by 
buoys  i  and  j . 
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The  p  (i/rr)  function^  if  the  two  recorded  sig- 
nals were  broadband^  of  long  duration  and  with  high  S/N 
ratioSf  would  appear  as  shown  in  Fig  U, 

Low  signal  to  noise  ratios  could  hide  the  real 
maximum  and  create  false  ones;  narrowband  or  coherent  sig- 
nals would  present  ambiguities  with  the  creation  of  many 
c 1 ose  max  i  ma . 

By  considering  the  availability  of  eguipments 
and  technigues  to  obtain  the  proper  elimination  of  the 
difference  in  freguencv  shifts  in  the  recorded  signals  ^ 
this  study  assumed  the  possibility  of  having  time  delay 
measurements  available  to  provide  information  about  the  tar- 
get position  according  to  the  relationship 


T.  .  (k)  =  td  .  Ck)  -  d.  (k)]  /  c   +   v(k) 
13  1        J 


(a. 37) 


where  v(k)  must  account  for  the  errors  caused  by  the  dif- 
ferent noise  contributions  in  each  signal/  for  the  accuracy 
and  orecision  of  the  cross-correlation  algorithm  or  devicer 
and  for  variations  in  t a rget -sensors  geometry  during  the 
recording  of  the  signals. 

Inaccuracies  anywhere  from  0.5  sec  down  to  a  few 
milliseconds  seem  very  reasonable  to  expect  in  many  cases. 

d.   Signal  Strength  Difference 

For   the   same   reasons   already   explained    in 
llfDtlfC       the   information   contained   in   the  difference  of 
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i|»U.-) 


(time  delay) 


(compression  factor) 


Figure  4  :  Correlation  function, 
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strength  of  the  collected  signals^  while  of  great  use  for  an 
experienced  and  skilled  audio  operator/  was  not  incorporated 
in  this  study . 


ai 


III.   THEORETICAL  SOLUTION 


A.   RECURSIVE  BAYES  FORMULATION 


The  maximum  amount  of  information  at  time  t,  f  in  a   pro- 
le 

babilistic  senses  is  concentrated  in  the  conditional  density 

function  of  the  states  given  all  the  observations  up  to  this 

k  k 

timer  p(x_(k)jZ  )/  where  Z   represents  the  set  {^(k),  z(k-l)/ 

...,  2(0)>  . 


By  specifying  a  cost  function  one  determines  how  this 
information  is  used  to  obtain  the  best  estimate  [5].  For 
example/  the  minimum  variance  estimate  is  the  mean  of  the 
conditional  density. 

From  the  assumptions  made  in  Section  II/A,  the  random 
processes  x(k)  and  {x(k)»z(k)}  are  Markov  and  the  condi- 
tional density  can  be  written  in  recursive  form  using  Bayes' 
Lawf  as  is  shown  in  (5] r      f6) f      17] f      [8] . 


p(x(k)|Z^)  =  c(k).p(x(k)  iz'^"^).p(z(k)  lx(k)) 


(3.1) 


where 


p  (  X ( k ) ! l^'h 


p(x_(k-l)  !Z^"1)  .p(x_(k)  !^(k-l)).d^(k-l) 

(3.1a) 


ana 


a2 


1  /  c(k)  =  p(z(k) IZ^"^  = 


p(x(k)  !Z^-^.p(z(k) !x(k)).dx(k) 


(3.1b) 


Given  the  initial  condition  p(x(0)|Z~-^)  =  p(x(0))  and  the 
probability  law  of  the  system  random  input  and  measurement 
noise#  p(w^(k))  and  p(v^(k))rone  can  proceed  to  evaluate  these 
equations  and  the  filtering  problem  can  be  regarded  as  hav- 
i  ng  been  solved  when  p(x^(k)iZ  )  can  be  determined  for  all  k. 

For  any  situation  different  from  the  one  where  the  sys- 
tem and  measurement  equations  are  linear  and  the  apriori 
distributions  are  Gaussian,  the  density  p(2<.(k)!Z  )  is  not 
Gaussian  and  generally  cannot  be  determined  in  closed  form. 


In  the  problem  studied  here^  the  conditional  densities 
in  Equation  (3,1)  are  nonlinear  functions  of  the  states  and 
the  measurements.  Linearizations/  Taylor  expansions  and 
Gaussian  aoproximations  will  be  used#  however,  as  an 
engineering  comoromise,  in  order  to  avoid  complex  numerical 
integrations  in  n  dimensions.  This  aoproach  was  chosen  tak- 
ing into  account  the  limitations  on  processing  eauioment 
available  for  the  practical  solution  of  the  addressed  track- 
i  ng  prob 1  em. 


43 


B.   ONE-STEP  PREDICTION 

Let  us  assume  that  at  time  t^  /  during  a  filtering  pro- 
cesSf  an  approximately  Gaussian  density  p  (  x  (  k- 1  )  !  Z '^"l)  has 
been  obtained  with  mean  value  p  (k-1)  and  covariance  matrix 
M(k-l). 

If  one  chooses  the  mean  of  this  conditional   density   as 


the  estimate  at  time 


V] 


x_(k-l)  =  E  [x.f  k-1)  ;Z^"^  =  y(k-l) 

the  estimation  error  at  time  k-1  has   the   same   probability 

k-L 
density  as  p(x(k-l)|Z  1    but  with  zero  meanr  i.e.f 


E[e(k-l)l  =  El(i{k-1)  -  x(k-l))|Z^"^  =  0 


(3.2) 


ECe(k-l)eMk-l)]  =  P(k-l)  =  M(k-l) 


(3.3) 


At  time  t,  a  new  set  of  measurements  arrives  and  a  pred- 
k 

iction   at   that   time   must   first  be  obtained  based  on  the 

measurements  up  to  t 

k-1 

The  dynamics  of  the  system  are  described  by  Equation 
(2.3)  and  it  can  be  seen  that/  even  if  p  (  x  (  k- 1  )  I  Z*^"-^)  were 
really  Gaussian^  the  density  p(x^(k)iZ  )  would  not  normally 
be  because  of  the  non  1  i near i t i es  involved.  Nevertheless  it 
will  be  assumed  that  a  reasonably  close  approximation  to  a 
Gaussian   density   exists   and   proceed   to  find  approximate 

values  for  the  first  two  moments  of  the  conditional  one-step 

k- 1 
pre-diction  density  p(x(k)|Z   ). 


aa 


The  first  step  is  to  exoand  Equation  (2.3)  in   a   Taylor 

series  about  the  estimate  xCk-l).  The  dependence  on  t,   and  t 

—  k       k-1 

k-1 
is  dropped  for  simplicity.   Note  that  Z    i^  known   and   was 

used  in  determining  x^(k-i). 

x(k)  =  f(x(k-l))  f  g(x(k-l)).w(k-n  = 


=  f(i(k-n)  +  |£ 

—  —  dX 


. tx(k-l )  -  x(k-l)J  +  ...  + 
i(k-l) 


+  a(x(k-l)) .w(k-l )+...= 


=  j^(£(k-n)  -  I 


n   n 
i=l  j=l 


n 

1 

i=l 

3f 

8f 

3x. 

1 

.e    (k-1) 

i 

x(k-l) 

t 
r  1^  -  1  "» 

3x. Sx. 
1      J 

•  c.vis       I.  /   »  <::       VIS       i/ 

x(k-l) 

+  . . .  + 


+  q(i(k-l  ))  .w(k-l  )  - 


n 

I 
i=l 


3  8. 


3x 


.w  (k-1) 
i 

x(k-l) 


.e( k-1 )  +  . . 


(3. a) 


where  x.  and  e.  are  individual  components  of  _x  and  (_x  -   x); 

g.    is  a  column  vector  of  g  /  and  w   is  a  component  of  w. 
1  —        i  — 

1.   Linear  Approximation 

In  this  case  only  the  terms  in  the   expansion   which 
are    linear  in  e(k-l)  and  w(k-l)  are    considered^  i.e.r 

x(k)  =  f(i(k-l))  -  ({)(k-1)  .e(k-l  )  + 


+  q(x(k-l)) .w(k-l) 


(3.5) 


as 


where 


(j)  (k-1  )  = 


_  9 


dx k   k-1 


(3.6) 


From  Equation  (3.5)  the  mean  value  of   p(x(k)IZ  ^  ), 
which  is  chosen  as  our  predicted  value/  is 


(k|k-l)  =  E(x(k)|Z^^  =  f(x(k-l))  - 


-  (|)  (k-l).E  te(k-l  )1  +  g(J(k-l))  .E  [w(k-nj 


But  e(k-l)  and  w(k-l)  have  zero  means^  therefore^ 


i(k|k-l )  = 


f  (x(k-l),t,  ,t,  J 

—  —        k   k-1 


(3.7) 


Thus  the  one-step  prediction  error  has  the  mean 


E[e(k;k-1)]  =  Et(;(k!k-1)  -  x(k));^~lj  = 


=  E  l(f,(k-l)  .e(k-l  )  -  g(x(k-l))  .w(k-l)J  =  0 


The  prediction  error  covariance  matrix^  which  is 
also  the  covariance  matrix  of  p(x(k)|Z  )/  is^  using  simpli- 
f  i  ed  not  at  i  on  » 


P(k;k-1)  =  E  [e(k|k-l)e'-(k|k-l)]  = 


C 


=  E  [(()(k-t).e(k-l)  -  g(x(k-n  )  .w(k-l)]  ( 


=^1  = 

.  .  •  J  I  ~ 


T   ,  T 


T  _T 


T  .  T 


.T  ^T 


=  E  (  ((),  e_.  e_  .  (j)  -  (j).e_._w  ._q   -  g.w.e  •({)   "♦"  g.w.w'-.g-'-l 


ab 


From  the   assumptions   about   the   independence   and 
discrete  white  noise  characteristics  of  Wf  the  terms 


((,  (k-1  ).Ere(k-n  .w^^Ck-l)]  .g'^(J(k-n)  =  0 


q(i(k-l)).E  tw(k-l).e'^(k-l)l  .4)'^(k-l)  =  0 


and 


P(k{k-n  =  (})  (k-l).P(k-l).  <j)Mk-l)  + 


+  q(x(k-l)).Q(k-l).gMx(k-l)) 


(3.8) 


2 .   First  Order  Approximation 

In  this  development  one  retains  all  terms   in   Equa- 
tion (3.4)  up  to  first  order  partial  derivatives^  i.e./ 

^(k)  =  f(x^(k-n)  -  (j)(k-l)  .e(k-l)  + 

+  a(x(k-l))  .w(k-l  )  - 


m 


y   A  (k-n  .w  (k-l)  .  e(k-l) 


i=l 


whe  re 


(3.9) 


A.  (k-l)  =  -^   a.  (i(k-l),t,  ,t   ) 
1         dx   — 1  —        k   k-l 


(3.10) 


As  in  the  last  section^  the  predicted  value  is 


^(klk-l)  =  f_(i_(k-l))  -  <^  (K-l)  .E  [e(k-l  ))  + 


m 


+  q(x(k-l) ) .E  (w(k-l )J  -  y   A  .E[w  ( k- 1 ) . e ( k- 1 ) ] 


lil 


From  the  assumptions  about  ^(k-1)   and   e(k-l)^   the 
predi  c t  i  on  is 

^(klk-l)  =  f^Cx^Ck-D), 

the  same  as  in  Equation  (3.7), 


The  one-step  prediction   error   has   zero   mean   and 
covariance  matrix  given  by/  in  simplified  notation/ 


[■ 


P(k|k-1)  =  E  [(|,(k-l).e(k-n  -  g(x(k-l))  .w(k-n  t 


m 


+   y   A.  (k-n  .w.(k-l).e(k-l)l  ( 

^     X         1        — 

r=i 


=  E  [  (j).  e_.  e^-"-  .  (J) 


T  ^T  . 


m 


•  •  •  J   I  ~ 


J   aT 


e.w-'-.q-'-  +   y   (b.e.e.A-'-.w 


m 


-g.w.e.(b   •♦•g.w.w'r.gT-   V    g.w.e^.A^.w   + 
"      i£l ii 


m 


T   T 


m 


.T  ^T 


+   y   A,.e.e  .A  .w.  -   y   A  .e.w  .q-'-.w   + 

111    ni  T     T 

+   y  y    A..e.e.  A-'-.v^.-wJ 


From  the  assumptions  about  w^  one  has  that 


t  (e(k-l) .w^(k-l )]  =  0 


E  (e(k-l)  .e  (k-1  )  .w.  (k-1  )1  =  E  [e  (  k-1  )  .eM  k- 1  )  1 

~       ~        1  —      — 


E  [w  (k-1)]  =  0 
i 


E  [w(k-l  )  .e  "-(k-l  )  .w  (k-1)]  =  E[w(k-n.w  (k-1)] 

i  —        i 


E  [e^(k-l )]  =  0 


as 


t[e(k-l).e(k-U.w.{k-l).w.  (k-1)]     = 
-  -  1  J 


=    E  [e(k-l)  .e    (k-1)]  .E  Iw.  (k-1 )  .w.  (k-1)]     = 
-  -  1  J 


=    P(k-l)     .    q. . (k-1) 
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where    q,,(k-l)     is    the    j^lth     term    of    Q(k-l) 


Thus  f 


P(klk-l)    =  (f)  (k-l).P(k-l).(J)  ^  (k-1)     t 


+    £(_x(k-l)).Q(k-l).jMx^(k-l))    f 


m     m 


I     I        f^    .P(k-l)  .A.  ,a.  .   (k-1) 


=lj=l 


J        Ji 


(3.11) 


If  the  Q  matrix  is  diagonal*  a.   =0  for  j  t    1  and 


P(klk-l)  =  (t)(k-l)  .P(k-l)  .,|,Mk-l)  f 


+  g(x(k-l)  )  .Q(k-l)  .aMx(k-l  )  )  + 


m 


I  A..P(k-l).ASa..(k-l) 


i=l 


T 

i   11 


(3.12) 


3 .   Other  Approximations 

Depending  on  how  many  terms  of  Equation  (3.^)  are 
considered/  other  different  prediction  values  and  prediction 
error  covariance  matrices  are  obtained.  Higher  moments  of 
the  conditional  densities  would  be  required  with  increased 
complexity  and  computational  burden. 
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C.   UP-OATING  THE  ESTIMATE 

Let  us  assume  that  we   have   an   approximately   Gaussian 

k-1 
density   p(x.(tc)IZ    )   whose   mean   value   is  the  prediction 

At 

x_(k|k-l)  and  whose  second  central  moment  is   the   prediction 
error  covariance  matrix  P(k!k-1). 

A  new  set  of  measurements  is  given^  consisting  of  non- 
linear noisy  observations  of  the  statesf  represented  by  the 
rel at  i  onsh  i  p 


2_(k)  =  h^(j<(k),t  )  +  v_(k) 


It  is  now  necessary  to  process  z(k)  and  extract  the 
information  it  brings.  A  new  approximate  Gaussian  density 
p(x^(k)iZ  )  is  assumed  to  be  generated  whose  moments  are  the 
new  estimate  i^(k)  and  the  estimation  error  covariance  matrix 

P(KJ. 


This  new  density  is  given  by  Equation  (3.1)  or  by 


p(x_(k)!z'^)  =  p(£(k)  Ix^Ck)  )  .p(j^(k)  ;Z^"1)  /  p(z_(k)  iz^-i) 

(3.13) 


k-1 
p  (_x  (  k  )  I  Z   j  is  assumed  to  be  approximately  Gaussian  with 

moments  ^(k|k-l)  and  P(kik-n. 


Expanding  the  measurement  eauation  around   x(k{k-l)   one 


has 


Z(k3  =  h(i(klk-l))  +  1= 


.[x_ik)     -  ;i(k  I  k-1)]  ♦  ...  +  _v(k) 
x(k|k-l) 
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Retaining  only  the  linear  terms,    one  has 

z^(k)  =  h_U(k!k-l)  +  H(k).[x^(k)  -  ^(k!k-l)]  +  _v_(k) 

(3.14) 


where 


H(k)  = 


^  h(x(k|k-l),t,  ) 

3x  —  —  k 


(3.15) 


From  Equation  (3.14)  one  obtains 


E  I  z  (  k )  I  Z^'h     =  h(i(k|k-l)) 


and  it  is  easily  shown  that 


Var  (z(k)  12^  -h  =  H  (  k  )  .P  (  k  I  k-1  )  .hF  ( ic )  +  R(k) 


Also  from  Equation  (3,14), 

E[z_(k)  :x_f  k)l  =  Jh(^(k|k-1)  )  f  H(k).[_x(k)  -  £(k|k-l)] 

Var  [zik) !^(k)]  =  R(k) 

If  one  now  assumes  that  p  ( z_(  k )  1 2r~-^)  and  p(2(k)|x(k))  are 
approximately  Gaussian  densities  with  the  moments  given 
above*  then  Equation  (3,13)  becomes 


.7K 


1/2 


p(x^(k)!Z'^)  =  |H(k)  ,P(k!k-l)  .HMk)  +  R(k)|  / 

r./o  1/2  1/2 

/     ((27rf^^,  iP(klk-l)  (     ,jR(k)|        ]     .    exp[B] 


where 
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exp  [BJ     =    exp 


{■^[[' 


(x(k)    -    i(k|k-l)l^.p-'-(k!k-l).l...] 


[■- 


+     |[2(k)    -  ]^(£(k|k-l))    -    H(k).[x(k)    -    x(k!k-l)]]T, 

R" 


■E 


■■"■(k).  t...]        - 


fzCk)    -    h(J(k|k-l))J  ^. 


,-l 


[H(k).P(k|k-l).HMk)     +    R(k)l  -".  {...] 


(3.16) 


Simplifying    the    notationr     this       exponent       can       be       rear' 
ranged    in    the    foHowinq    way 

(x    -    i)'^,P~\(x    -    5)     +     tz    -    hCx)     -    H.(x    -    J)]T  .R-1.  [...]     - 


-     Iz    -    h()?)]'^.  [H.P.h'^+    Rr\[z    -    hC^)]     = 


=    (x    -    x)^P-'(x    -    i)     +    (x    -    i)^H^R"^(x    -    x)    - 


-     Iz    -    h(x)F    R"1h(x     -    x)     -     (x     -    xf  H^R^Cz    -    h(i)]     + 


■»•     [z    -    h(i)]^R^(z    -    h(S)]     - 


-     tz    -    h(x)]^   (HPH^     +    Rf-^lz    -    h(x)]     = 


-    T        -1  T   -1 

(X    -    x)     (P       f    H    R  -Vl)  (x     -    x) 


-      (  X      -      X 


T     T    -L 
y H^R   Iz    -    h(x)]     - 


-     tz    -    h(x)]^R^(x     -    x)     - 


[z    -    h(i)]^   [R    ■^-     (HPH^     +    R)^]   [z     -    h(x)] 
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Using  the  matrix  inversion  lemma  [91, 


R  ■•■  -  (HPh"^  tR)"^  =  R  H(P  ^+  H^ff^y^H'^R"^ 


Since  R  and  P  are  symmetric  matrices 


r-^CP"-'-  +    h'^R  ■•■H)  W  R    ■•■  =    A-^lP  ■"■   +    hf  R'-'hI  a 


where 


A    =     (P  ■"■  +    h'^R-'-Hf-Si^R"-'- 


and  the  exponent  of  Eauation  (5.16)  becomes 


-  T   -1    T  -1 
(x  -  x)(P+  H^R-^H)(x  -  x)  - 


(  X  -  X  )^  H^  R  ■"■[  2  -  h  (  i  )  1 


-  [2  -  h(x)l  R  \Mx    -  x)  - 


-  [(P"^+  H^R^hT-'-H^K  ^(z  -  h(i))j^lP"^+  H^ffVl][...) 


which  is  a  quadratic  form  that  can  be  expressed  as 


X  -  X 


-  [P"-^  +  H^  R'h]  \^^Rhz    -  h(i) 


.  ( P  1  +  h'^  R"^H  ]  . 


V-   _J 


And  finally  the  conditional  density  of   Equation   (3.16) 
can  be  put  in  the  form 


D(x(k) !Z^)  =  a(k) .exD 


-|-[x.(l<)  '    i_ikn\p~\<)  .[.,,] 


where,  by  definition, 
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x(k)  =  iCklk-l)  +  P(k)  .HT(k)  .R-^k)  .  [z(k)  - 


h(x(k|k-n,t  )] 

—  —  k 


(3.17) 


and 


P"\  k  J  =  P"\  k  I  k  -  1  )  +  H^  (  k  )  R  \  k  )  H  (  k  ) 


or#  using  the  matrix  inversion  lemma  again/ 


P(k)  =  P(klk-l)  -  PCklk-DH^Ck).  (H(k)P(k|k-l)H^  (k)  f 


+  R(k)r-'-.H(k)P(k|k-l) 


or 


P(k)  =  (I  -  G(k)H(k)] .P(k|k-1 ) 


(3.18) 


where 


G(k)  =  P(k;k-n .H^ (k) . rH(k) .P(kik-i) .HT(k)  +  R(k)] 

(3.19) 

From  hquation  (3.19)  one  obtains  - 
G(k)  =  P(k!k-1) .H^ (k)  .R"^k)  . 

[H(k).p(kik-i)  .H'^(k).R'\k)  +  ir^ 

Premu 1 t i p 1 y i ng  by  P(k)P  (k)/  and  using  the  expression   found 
earlier  for  P  vk)  gives 


G(k)  =  P(k).[I  +  H^R"-'-HP(k|k-l)]  H^R  \ 


(HP(klk-l)H^R  ■'■  +  I]  ^  = 


5a 


=  P(k)  .  (H^R  ■"■  +  H'^R-^HP(k|k-l)H^R  -^  . 

(I  +  HP(k;  k-nn'^R"  ■')"■'■ 


and  final  1 y » 


G(k)  =  P(k).H'^(k).R  \k) 


(3.20) 


Applying  this  result  to  Eauation  (3.17)  gives 


i(k)  =  i(klk-l)  •»-  G(k).(2(k)  -  h.(_x(klk-l),t^)] 


(3.21) 


If  more  terms  were  taken  in  the  expansion  of  h(x(k)#t,  ) 

—  —      k 

than  the  ones  considered  in  Equation  (3.1^)»  different  esti- 
mates and  covariance  matrices  would  be  obtained.  Higher 
moments  of  the  conditional  densities  would  then  be  reauired. 

We  are  now  in  position  to  restart  the  one-step  predic- 
tion of  Section  III^B  and  process  a  measurement  ocurrinq  at 
time  t^^j^. 

The  initial  step  in  the  whole  process  is  at  t   where   it 

o 

is  assumed  that  an  approximate  Gaussian  distribution  exists 
for  the  initial  value  of  the  statesr  ^(O),  Choosing  the  ini- 
tial  estimate   as   the  mean  value  of  this  distribution^  one 

has 

x_(  0  )  =  E  [j(  (  0  ) ) 

E  le(0)l  =  E  [i(0)  -  x(0)]  =  0 
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P(0)  =  ECe(O)eMO)]  =  Varix^(O)] 

The  set  of  Equations  (3.7),  (3.8),  (3.18),  (3.19)  and 
(3.<fl)  is  generally  known  as  the  Extended  Kalmam  Filter 
equat  i  ons . 

D.   PARAMETERS  FOR  THE  TRACKING  PROBLEM 

In  order  to  apply  the  above  relations  to  the  passive 
tracking  problem,  a  series  of  vectors  and  matrices  first 
must  be  determined. 

1  .   System  Dynami  cs 

Applying  definition  (3.6)  to  Equation  (2.10)  yields 


(()(k)  = 


(3.22) 


At .cos  b ( k ) 

At .sin  b(k) 

1 


-  At.s(k).sin  b(k) 


At .s(k) .cos  b(k) 


where   At  =  ( t,  ,  -t,  ) 
k+1   k 


Breaking  up  Equation  (2,11)  into  its  column   vectors 


gi  ves 
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g  (x  Ck)  , t   ft  )  = 

-1  -     k+1  k 


1  At^  .cos  b(k) 

2 

iAt^  .sin  b(k) 
2 


At 


(3.^3) 


g  (  X  (  k  )  f  t,   r  t ,  )  = 

-2  -  •     k+1   k 


l-At^  .s(k)  .sin  b(k) 
2 


^At^ .s(k).cos  b(k) 
2 


At 


(3.2a) 


Definition  (3,10)  requires 


A  (k) 

1 


(3.25) 


0    0    0 


0    0    0 


0    0    0 


0    0    0 


^At^  .sin  b(k) 


lAt^  .cos  b(k) 
2 
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A2(k)  = 


(3.26) 


0    0 


0    0 


0    0 


0    0 


-  ^  A  t^  .  s  i  n  b  (  k  ) 


^ A  t  .cos  b(k) 


^At^  .s(k) .cos  b(k) 


^At^  .s(k).sin  b(k) 
2 


Assuming  independent  random  forcing  inputs  yields 


Q(k)  =  E  [w(k)w'-  (k)]  = 


a^(k)    0 
s 


0   c^Ck) 
b 


(3.27) 


Equations  (3.8)  and  (3.11)  require  the  product 


g(x(k),t   ft  ) .Q(k) .gT(i(k) ,t   ,t  )  = 
k+1  k        -  -      k+1  k 


=  At^  . 


a  (k)     0 
s 


aj(k) 

D 


(3.28) 


where 


a  =  7  At^  .  (cos^  b(k)  .  a^(k)  +  s^  (  k  )  .  s  i  n^b  (  k  )  .^  ^^(  k  )  ] 

T'  S  D 

b  =  ^At^.sin  b(k).cos  b(k).[a^(k)  -  s^  (  k  )  .^  ^(  k  )  ] 
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c  =   7At ^. (sin^b(k) . a^Ck)  +  s ^ ( k ) .cos^b ( k ) . a^ ( k)  ] 

^  S  D 


and 


1  "2 

d  =   -rAt.COS  b(k).a  (k) 

2  s 


1     -^  '^2 

e  =  --rAt.s(k),sin  b(k),a,  (k) 
^  b 


1  "2 

i  =   ^At.sin  b(k),a  (k) 

2  s 


j  =   ^At.s(k).cos  b(k).a^(k) 
2  b 


If  the  expanded  state  variable  formulation  is   usedf 
the  matrices  given  below  are  needed 


Q®(k)  =  Elw^(k)w®^(k)]  = 


Q(k)     !  0 


0 


0     0 


■-I— 
!   2 


(3.29) 


g  .Q  .g    = 


g.Q.g 


,2   2 

At  •  Of 


(3.30) 


*^  = 


0     '  1 


(3.51) 
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2.   Measurement  s 

In  the  following  development  Equation  (3.15)  is 
applied  to  each  of  the  measurements  described  in  Section 
IlfD  to  yield  the  appropriate  linearized  observation  matrix. 

a.   Bearing  Measurement 

From  Equation  (2.2^)  one  has 

8  b.  =  (x  -  x.)^/[(x  -  X.  :?+  (y  -  y.)^]  . 


55 


^  [(y  -  y.  )/(x  -  X.  J] 


and  then 


T.i 


(k)  = 


^b.  (xCk|k-l),t,  ) 
dX   X  —  k 


=[- 


[y(k;k-l)  -  y.(k)]/a      [x(k|k-l)  -  x 

1 


.  (k)]/a 


0    0 
(3.32) 


where 


a=  d^(k)  =  (x(klk-l)  -  X  (k)]^   ^-  [y(klk-l)  -  y  (k)]^ 
1  i  i 


(3.33) 


and  d.(k)  is  the  predicted  distance  of  the  target  from   buoy 


b.   Frequency  Measurement 


From  Equation  (2.25)^ 
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3C      1 


=     (c/(c    - 


s.)).Af       +     [f    .c/(c    -    s.)^].-i.s 


ac    o 


3C      1 


f   /f 

i       o 


8?      o 


f^  /(f      .c)      .   -iL-   S 


3?      i 


where    s       is    given    in    Equation    (2.26)    and 
i 


3  3 

V- s  .    =    -    cosCb    -    b.).-7—  s    + 
8?      1  1        8C 


+    s.sinCb    -    b    ).-— rb    " 

1        3C 


-    S.sinCb    -    b    ).-— rb 

1  dC       1 


and    t  henf 


(3.3^) 


H^.  (k) 

fi 


-;7-f.    (x(k|k-l),t,    ) 
3x     1     —  k 


=     [y(klk-l)-y.(k)]/c         - [5 ( k I k- I ) -x . ( k ) ] / ^       6     y 


.'••] 


where 


(3.35) 


f     (x(klk-l)  ,t ,) 
1    —  k 


f     (klk-1)  .c    /     [c    + 
o 


+    s(k| k-1) .cos(b(k ! k-1)    - 


b    (k))J 

i 


b . ( k )    =    arc  tan 
1 


[y(k;k-l)    -    y     (k)l  / 

1 


[x(k! k-1 )    -    X 


J.))] 


(3.36) 


Y=    f^  .s(k!k-l)  .sin(b(k;k-l )    -    b    (k))    / 
i  i 
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(f   (kik-n.c) 

o 


3  =  -  f  %cos(b(k!  k-1)  -  b  (k))  /  (f  (klk-n.cJ 
1  i  o 


6  =  a/ Y   '  a  as  in  (3. 33) , 


c.   Frequency  Difference  Measurement 


From  Equation  (2,2b)f 


-^  (f.  -  f  .)  =  f  /c  .  -1.  (s  .  -  s.  ) 
3C   1     J      o      3C    1     J 


where  —  s   is  given  in  Equation  (3.3^).  Then^ 
3C   i 


H  ^.(k)  =  J-  (f   -  f  ) 
di- 


3x    i 


=  f  /c  .  (h    h    h    h  ]  (3.37) 

o        12    3    4 
x(k|k-l) 


where 


h^  =  a  .[y(klk-l)  -  y  (k)]  -  a  .[y(k|k-l)  -  y  (k)) 

h^=  a  .tx(k!k-l)  -  X  (k)3  -  a  .[x(k|k-l)  -  x  (k)] 

2  j  j         i  i 

h^  =  cos(b(k|k-l)  -  b  (k))  -  cos(b(k|k-l)  -  b  (k)) 

3  j  i 

h,  =  s(k |k-l)  .  [sin(b(k!k-l )  -  b  (k))  - 

4  i 


sin(b(k!k-l)  -  b  (k) )] 


a   =  s(k| k-1) .sin(b(k ! k-1)  -  b  (k))  /^ 
■  1  i       "^ 


b.(k)  as  in  (3.36)  and  ct  as  in  (3.33). 
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d.   Frequency  Ratio  Measurement 


From  Equation  (2.28)^ 


-L  (t.  /f  )  ;  i/c  .  J-  (s  -  s  ) 


yC 


n     i      j 


Comparing  with  the  frequency  difference  measure 
mentf  one  has 


H   (k)  =  -£-  (f  /f  ) 

ri         32i    i   3 


=  1/f   .  H   (k) 

o     di 

x(k|k-l) 


(3.38) 


e«   Time  Delay  Measurement 


From  Equation  (2.36)/ 


T       =    1/c     .    ^(d       -    d    ) 


2c      ij 


3?        i  j 


where 


d2  =    (x    -    X    )2     t    (y    -    y    )2 


and    then/ 


h^,.(k)    = 


-At.. 
ax    ij 


=    i/c    .    Ih        h  2     0        01 


x(k|k-l) 


where 


h  ,   =     [x(l<  |k-l)    -    X  .(k)l     /    d  .(k)    - 
1  11 


(3.39) 


tx(klk-l)     -    X  .(k)]     /    d  .(k) 
J  3 
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h   =  ly(k|k-l  )  -  y  (k)]  /  d  (k)  - 
^  i  i 


[y(klk-l)  -  VjCk)]  /  d .(k) 


and  d   is  given  in  Equation  (3,33), 
i 


6a 


IV.   SPECIAL  CHARACTERISTICS  OF  KALMAN  FILTERS 


Equations  (3,7)  and  (3.8)  characterize  the  one-step 
prediction  of  an  Extended  Kalman  Filter;  Equations  (3.18)» 
(3.19)  and  (3.21)  represent  the  estimation  update. 

iwhen  implemented  in  a  computing  device^  the  timing  of 
these  steps  is  shown  in  Fig.  5.  The  actual  state  is 
represented  by  a  broken  line  and  the  output  of  the  filter  by 
the  solid  line.  In  general*  one  is  dealing  with  vector 
processes  so  Figure  5  is  representative  of  only  one  com- 
ponent. Nevert hel  ess #  this  reoresent at i on  is  helpful  in 
visualizing  the  steps  involved  in  the  estimation  process. 

In  a  general  problem  the  times  of  occurrence  of  meas- 
urements  are   not   eaually   separated   and  are    not  known  in 

advance.  At  a  time  t  r    when  a  new   set   of   a   measurements 

k-1 

becomes   available*   the   output   of  the  filter  is  still  the 

last  estimate  made  close  to  t   *  unless  predictions  are    made 

k-2 

at  regular  time  intervals  independent  of  the  spacing  between 
measurement  s . 


A  time  period  T   is  spent   in 

P 


a  prediction  phase  to 
evaluate  Equations  (3.7)  and  (3.8)*  after  which  the  best 
information  about  the  state  of  the  plant  is  the  oredicted 
value  x(k-l!k-2). 
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If  the  times  of  occurrence  of  measurements  were  known 
in  advance^  the  predictions  could  be  made  any  time  the  sys- 
tem was  idle  waitinq  for  the  new  measurements  and   T    could 

P 

be  made  zero.  ' 

Nextf  the  q  measurements  are  processed  using   Equations 

(3.13)r   (3,19)   and   (3.21)  and  a  time  T   is  spent  which  is 

e 

relatively  large  mainly  because  of   a   matrix   inversion   of 

dimension   (q  x  q)f  in  the  gain  equation^  (3.19).  Only  after 

this  time  has  passed  is  the  new  estimate   x(k-l)   available. 

This  estimate  will  be  the  output  of  the  filter  until  t,   when 

k 

the  process  is  repeated. 

Some  important  points  in  this  estimation   process   must 

be  stressed: 

---  the   new   estimates   are   only   available   T    +   T 

P       e 

seconds  after  the  arrival  of  new  measurements. 

---  a  major  oart  of  the  time  oeriod  of  duration   T    is 

e 

spent  in  the  gain  equation  due  to  the  required  matrix  inver- 
sion. 

---  from  the  ena  of  the  computation  of  an  estimate 
until  the  arrival  of  a  new  set  of  measurements  the  computing 
equipment  is  idle  and  thus  free  to  execute  other  chores. 

---  the  one-step  prediction  is  based  on  a  linearization 
of   the  plant  dynamics  over  an  extended  periods  for  example* 

from  t    to  t  . 

k-1     k 

---  the  measurement  equations  are  based  on  a  lineariza- 
tion about  the  predicted  values. 
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A,   PARTITIONING  OF  MEASUREMENTS 

When  some  of  the  measurements  come  from  statistically 
independent  sources/  steps  can  be  taken  to  alleviate  orob- 
lems  which  arise.  Considering  the  special  case  when  all  the 
measurement  noise  components  are  independent*  the  covariance 
matrix  R(k)  is  diagonal  and  the  following  development 
app lies. 

1  •    The  L  i  near  Case 


Assume  a  dynamic  system  with  linear  observations   of 


the  form 


2(k)  =  HCk) .x(k)  +  v(k) 


(a.i) 


where  HCk)  is  not  a  function  of  the  states. 

The  estimation  equations  of  Chapter  III   become   the 
standard  Kalman  Filter  equations 

G(k)  =  P(k!k-l)H'^(k)  [H(k)P(k|k-l  )HT(k)  +  R  (  k  )  F^      (a. 2) 


J(k)  =  iCklk-l)  +  G(k)[2(k)  -  H(k).x(k;k-l)l 


(a. 3) 


P(k)  =  II  -  G(k)H(k)]P(k|k-l) 


(^.a) 


The  measurements  are  assumed  to  occur  simultaneously 
andf  as  shown  for  the  first  time  in  19],  for  this  linear 
case,  the  results  obtained  by  processing  all  q   measurements 
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at   once   using   these   eauations   are  the  same  as  the  final 
results  obtained  by  processing  one  measurement  at  a  time. 

In  the  latter  approach^  the  result  of  processing  one 
measurement  component  is  used  in  the  following  computation 
to  process  the  next  measurement  component. 

A  general  way  to  show  that  the  two  approaches  yield 
the  same  estimate  after  all  measurement  components  have  been 
processed  is  given  in  the  following  paragraphs. 

In  short  it  is  to  be  proved  that  the  estimates  x(k) 
and  P  ( k )  obtained  by  the  use  of  Equations  (4,2)^  (4.3)  and 
(4.U),  when  all  measurements  are  grouped  in  a  vector  z^,  are 
the  same  as  those  obtained  by  the  use  of  the  following 
iterative  equations  q  times/  once  for  each  component: 


G  .    =    P.    ,.H  ."^Z     fH.  P.      H^    f    R.  ) 

X  1-1  1  1     1-1       1  1 


(a. 5) 


X.     =    X.        +    G.  [z.     -    H.  X.      ]     = 
— 1         —1-1  1  —1  1—1-1 


=     [I-GH]x  +Gz 

i    i    —  i-1  i   i 


(4.6) 


P.      =     [I     -    G.  H  .1  .P        ■ 
1  11        1-1 


(4.7) 


i     =    1,     2#     .../     q 


where 


P       =    P(k|l<-1), 


X       =    x(klk-l) 

— o         — 
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P   =  P(k), 

q 


X   =  x(k) 
-q     - 


z(k)  = 


"^  2      " 

1 

r  H  ~ 
1 

'2 

H2 

• 

= 

• 

• 

• 

2 

H 

_   q_ 

_  q- 

x(k)  ♦ 


V 

L  q 


and   E[v  1  =  R^ 
i 


If  Equation  (^.6)  is  applied  a  times*  the   following 

relationship   between   x(k)  (or   x   )  and  x(k|k-l)  (or  x  ) 

-  -q        -  -o 

pesul t s 


x(k)  =  [I  -  G  H  ]  [I  -  G   H   ]...ri  -  G  H  ]x(k|k-l)  + 
-  q   q         q-1  q-1  11- 


+Gz        +CI-GHIG7  +...+ 

q      q  q    q       q-1   q-1 


+     [1    -    GqH^tl    -    Gq_3^H       ]...[!    -    G^H^JG,  z 


2    2'     I'l 


(a. 8) 


and  recursive  use  of  Equation  (^.7)  gives 


P(k)  =  [I  -  G^H^l  II  -  G  , H   J...CI  - 

q  q      q-1  q-i 


G^H^]P(k|k-l) 


(a. 9) 


Comparing  (a, 8)  to  (a, 3)  and  (a. 9)  to  (a. a)  one  sees 
that  the  assertion  can  be  proved  bv  showing  the  validity  of 
the  re  1  at  i  ons 
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II    -    G(k)H(k)l     = 


II    -G    H    ]  (I 

q    q 


G       H       1  ...  (I    -    G    H    1 
q-1    q-1  1    1 


(a. 10) 


G(k)    = 


(n    X    1) 


II    -    G    H    ]...[!    -    G^H^JG, 
q    q  2    2      1 


(n    X     1) 


[n    X     1) 


(I    -    G    H    IG     , 

q   q     q-1 


(^.11) 


OTf     in    3    compressed    way* 


G(k)    = 


(n    X     (q-t)) 
[  I    -    G    H    1  G* 

q   q 


(n    X     1) 
G 


(a. 12) 


and 


[I    -    G(k)H(k)]     =     [I    -    G    H    1[I    -    G*H*1  (a. 13) 

q  q 


where  G  is  the  gain  obtained  from  Equation  (3.19)  if  one 
has  a  q-1  dimensional  measurement  vector/  and  H  and  R  are 
equivalent  matrices^  as  shown  below 


H(k)  = 


(a.ia) 


R(k)  = 


(a. 15) 
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,*T 


*T    .   n,*,-l 


G*  =  PClc!k-l)H""   tH'^PCklk-l  )H''^  +  R"l 


(a. 16) 


Let  us  expand  Eauation  (^.2)  and  try   to   manipulate 
it  into  the  form  of  Ecuations  (4.10)  -  (4.13).  We  have 


[HPH   +  RJ  = 


Hj,PH^+  R^ 


H2PH^ 


T 
H  PH. 

q  1 


H^PH^ 


HjPH^*  R^ 


T 
H,  PH^ 

1  q 


H_PH 
2    q 


H  Phf  +  R 

q  q   (^ 


or 


[HPH   +  R]  = 


T 
^3 


^2 


(4.17) 


where 


A  ,  =  H*PH*T  >  R* 


i-2 


h*ph'^ 


^3 


H  PH*^  =  a^ 
q       —2 


H  PH^  +  R 

q  q    q 


Next  define  the  scalar  c  as 


c  =  a   -  a^  _  £ 
3  12 
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ThuSf 


c  = 


H  PH^  +  R   -  H  PH*^  [H*PH*'^+  R*r^H*PHT 
q    q      q      q  q 


Using  the  definition  of  G   in  Equation  (^,16)  gives 


C  =  H  PH^  +  R   -  H  G*H*PH^  = 
q    q      q      q        q 


=  H  (I  -  G*H*1PH^   +  R 

q  q    q 


(a. 18) 


It  is  shown  in  [10]  that 


\    1  12" 

T         i 
5-3        i        ^ 

-X 

B            !        b 
1         I       -2 

where 


«1  =  *!'  '    ''^-24^'="' 


^2  -  -  S^2^ 


Ca.l9a) 

(a. 19b) 


.T  _      T/-1  -1 
^3  -  "  ^3'l^ 


.   -1 


b,  =  c 


(a. 19c) 
(a. 19a) 


Using  these  rel at i onsh i psr 


G  =  PH*^  [HPh'^  +  R]  1  = 


=  tPH*'^   I  PH  J 


-^3 


^2 
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1  q-3 


PH*^b       +    PH^b    1 
-2  q     4 


=     [C     I     Dl 


(4.20) 


Applying    Equations    (a. 19a)    -    (a.l9d)    to    (a. 20), 
D    =    .    PH*TtH*PH*T     t    R*rlH*PHVlf    PRTc'lr 

q  q 

=    -    6*H*Ph'^c"'^+    PH^c'-'s     (I    -    G*H*1Ph'^c"^= 

q  q  q 

*     *         T  *     *         T  -1 

=     [I    -    G    H    IPH     [H     (I    -    G    H    )PH       +    R    1  -^ 
q       q  q  q 


Let  P*  =  II  -  G*  H*  JP  be  the  estimation  error 
covariance  matrix  after  the  processing  of  a-1  measurement 
components,  as  seen  from  Equation  (3.18).  Then, 


D  =  P*  H^CH  P*hT  +  R  r^  =  G 
q   q    q      q         q 


(a. 21) 


is  the  gain  when  processing  the  qth   measurement. 


'^e  also  have 


C  =  PH*^(H*PH*T  ^  R*rl  +  PH*T[H*PH*T   ^  r*]-1h*PhTh  PH*T 


q  q 


[H  PH  ^  +  R  rh'-^  '    PH^H  PH**^  (H*PH*'r  +  R*rVl  = 


q  q 


=  G*  f   g*h*phTh  gV^-  phTh  g*c-i 
q  q        q  q 


Since  c  is  a  scalar,  one  can  write 


C  =  G   -  II  - 


G*H*)PH^c"^H  G*  =  G*  -  DH  G* 


and,  therefore. 


7a 


C    =     [I    -    G    H    ] G 

q      q 


(a. 22) 


and    f  rom    C^  .  15) / 


.  =  [. 


G    =    I  tl    -    G    H    ]G 

q  q 


(a. 23) 


so 


It    is    also    the    case    that 


;h  =      I 


II    -    G    H    ]G  G 

q  q  q 


] 


H 


=     tl    -    G    H    1G*H*+    G    H 

q    q  q   q 


1-GH    =    I-GH-     tl-GH    ]G*H* 

q  q  q   q 


andf    thus 


I    -    GH    =     [I    -    G    H    ]  II    -    G*H*] 

q   q 


(4.24) 


Equations  (4.23)  and  (4.24)  are  the  same  as  Equa- 
tions (4,12)  and  (4,13)/  proving  the  initial  proposition. 

Computationally/  one  can  replace  a  (qxq)  matrix 
inversion  and  a  (nxq)  x  (qxq)  matrix  multiplication  by  only 
q  scalar  divisions.  For  an  indication  of  how  much  computing 
time  can  be  saved/  consider  that  a  matrix  multiplication 
requires  approximately  nq^  scalar  multiplications  ana  addi- 
tions/  and  that  a  matrix  inversion  requires  at  least  q  mul- 
tiplications and  additions/  not  counting  logic  and   indexing 
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time  nor  considering  storage  requirements. 

The    partitioning    of    the    measurements     then 

2 
corresponds   to   trading   at  least  q  (n  +  q)  multiplications 

and  additions  for  q  divisions.  For  a   machine   with   average 

multiplication   time   of  10  microseconds^  addition  time  of  2 

microseconds  and  division  time  of  12  mi c roseconds »  if  n   and 

q   are  equal  to  ^  a  saving  of  about  1,5  ms  is  obtained?  if  n 

and  q  are  of  the  order  of  20,  this  saving  could   be   of   the 

order  of  200  ms  at  each  measurement  time. 

Besides  this  comouting  time  saving,  the  processing 
of  one  measurement  at  a  time  provides  partial  estimates  that 
can  be  used  as  improvements  over  the  predicted  values,  as 
shown  in  Figure  6.  Also,  the  final  estimate  is  obtained 
sooner,  although  with  the  same  accuracy  as  when  the  measure- 
ments are  all  processed  simultaneously. 

2.   The  Nonlinear  Case 


We  assume  the  same  system  as  before   but   with   non- 
linear observations  of  the  form 

z(k)  =  il(x_(k),t^  )  +  v(k) 

By  linearizing  this  equation   as   shown   in   Chapter 
III,  using  the  definition 


H(k)  = 


_   3 


^  h(x(i<;k-i),t, )  , 

3x  —  —  k 


the  proof  of   the   last   section   is   still   valid   for   the 
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x(k-l) 


'x(k) 


x(k) 


x(k|k-l) 


\      \^^P 


measurements 
available . 


t  +T  +T 


time 


predicted  value 
available;  start 
of  processing  all 
q  measurements. 


P   e 


end  of  processing 
all  q  measurements 
estimated  value 
available . 


I    /x(k) 


—  -  —T 


x(k)  =  x(k) 

-^ 1 


x(k|k-l) 


A 

X(K- 


(k-:^) 


1 

\start 


predicted 
available ; 
of  processing 
first  measurement 


of 
processing 
second  meas 


time 

_end  of  processing  all  q 
measurements;  final  estimated 
value  available. 


Figure  6  :  Partitioning  of  measurements 
Linear  case. 
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linearized  equations. 

Iff  howeverr  the  partial  estimates  obtained  by  pro- 
cessing one  measurement  (hopefully  a  better  approximation  to 

x(k)  than  x(k|k-l)  results)  are   used   to   calculate   the   H. 

—       -   —  1 

vector  for  the  processing  of  the  n^xt  measurements  better 
accuracy  should  be  expected  than  when  all  H.  vectors  are 
obtained  using  the  relatively  crude  predicted  value.  This 
effect  should  be  particularly  significant  when  the  time 
between  each  group  of  measurements  is  relatively  large. 

Thus*  partitioning  the  measurements  in  an  Extended 
Kalman  Filter  should  provide  improvement  in  estimation  accu- 
racy in  addition  to  the  advantages  pointed  out  in  the  previ- 
ous section.  This  improvement  can  be  visualized  as  shown  in 
Fi  gure  7 . 

The  expected  improvement  in  the  estimation  process* 
brought  by  partitioning  the  measurements  in  Extended  Kalman 
FilterSf  can  be  explained  in  terms  of  the  conditional  densi- 
ties. 

Consider  two  independent  simultaneous  measurement s ^ 


z(k)  = 


z^  (k) 


Z2  (k) 


h  (^(k))  +  V  (k) 


h^  (x(k))  +  V  (k) 

2  -         2 


The  conditional  aensity   after   processing   the   two 
measurements  is^  from  Equation  (3.1)f 
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x(k) 


,J^' 


I  x(k) 


I 


x(k|K-l) 
I 

xCk-l) 


J— ^(k) 


\ 


estimate  if 
processing  the 
q  measurements 
together . 


^   \ 


+T 


end  |of 
processing  q 
measurements 
one  at  a  time 


time 


end  of  processing 
q  measurements 
alltogether . 


Figure  7  :  Partitioning  of  measurements 
Nonlinear  case. 
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p(j((k)|Z^)  =  pU(k)  !i  (k),2  (1<),Z^  ^  = 


=  p(z  (k),z  (k)  U(k)).p(j((k)  i/^^'^J/pCz  (k),z  (k)|Z^"-S  = 
J.      ^  X      ^ 


CpCz^Ck)  lz^(k),x(k))  /  p(22("<)lz  (  k  )  ,  z"^"^)  J  . 

(pCZj^  (k)  lx.(k)).pCj<(k)  !Z^"-^)/p(z^(k)  |Z^"^)1 


Butf  from  (3.1) 


pU(k)  ;z^(k),Z^  ^  =  p(z^  (k)  ;x(k)).p(x.(k)  iz'^  ^)  / 


p(z^(k)  !Z^  ^) 


al  sof 


p(Z2(l<)  !z^  (k)  ,x_(k))  =  pCz^  Ck)  |£(k)) 


because  of  the  i ndeoendence  assumption. 


And  sor 


(x(k)IZ^)  =  tp(z^(k)  IxCk))  /  o(z  (k)lz  (k),Z^  "'•)]. 
-  2     -  2       1 


pU(k)  lz^(k),Z^""S 


(a. 25) 


The  terms  on  the  right  hand  side  of  (^.25)  are: 

k-1 
p(x^(k)|z^  (k)/Z   )  ---  has  all  the  information  about  the 

states  from  processing  up  to  the  first  measurement  component 

z^(k)  . 


pCz^  (k)!£(k)) 


is   obtained   directly   from    the 
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measurement   equation   and  doesn't  depend  on  the  first  meas- 
urement . 

p( z  ( k ) I z  ( k )  / Z  )  ---  is  obtained  from  the  measurement 
equation  but  using  all  the  results  given  by  the  processing 
of  the  first  measurement  component. 

This  last  term  indicates  that  information  is  lost 
iff  in  a  Extended  Kalman  Filter^  the  linearization  for  pro- 
cessing a  measurement  is  based  on  the  predicted  value 
instead  of  on  the  value  estimated  after  processing  the 
preceding  measurement  component. 

3 .   The  Tracking  Problem 

In  the  tracking  problem  here  addressed^  the  measure- 
ments naturally  occur  at  different  times.  Bearing  indica- 
tions are  basically  continuous  or*  if  preprocessed r  succes- 
sive measurements  are  available  within  short  intervals  of 
time.  Frequency  indications  must  be  obtained  from  digital 
processing  of  records  of  considerable  time  length.  Time 
delay  indications  are  also  obtained  from  digital  processing 
but  with  longer  time  intervals. 

In  earlier  work,  til/  [21  and  [31/  these  measure- 
ments are  forced  to  occur  simultaneously  and  are  processed 
toget  her . 

In  the  previous  sections  of  this  chapter  it  was 
shown   thatf   if   the   measurements  are    all  available  at  the 
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same  timer  they  should  be  processed  seoarately  for  better 
computing  efficiency  and^  with  nonlinear  measurement s r  for 
possibly  better  estimation  accuracy. 

Since  one  is  now  free  from   having   to   process  the 

measurements   together/  and  encouraged  to  do  sor    why  not  use 

the  idle  time  of  the  computing  eguipment  and  process  the 
measurements  at  different  times? 

If  one  does  this  the  decrease  in  computation  time 
obtained  previously  will  be  diminished  because  of  the  extra 
computation  time  required  to  do  many  predictions^  one  for 
each  time  a  measurement  is  to  be  processed.  Nevertheless* 
the  advantages  are    ootentially  of  great  value*  that  isr 

---  a  measurement  can  be  processed  as  soon  as  it  is 
available^  with  no  waiting  for  other  measurements, 

---  measurements  that  occur  more  often  can  be  pro- 
cessed more  times  than  measurements  ocurring  less  fre- 
quent 1  y . 

---  the  output  of  the  filter  is  updated  frequently. 

---  the  one-step  prediction  is  based  on  lineariza- 
tions through  much  smaller  time-steps  than  before*  hopefully 
with  improved  accuracy  and  fewer  divergence  problems. 


Figure  8  gives  an  idea  of  the  improvement  one 
expects  to  obtain  by  (1)  processing  simultaneous  measure- 
ments separately  and  (2)  processing  the  measurements  as  they 
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naturally   occupf   when   a   Extended   Kalman  Filter  is  being 
used. 


B.   GRAPHICAL  INTERPRETATION 

General  understanding  of  the  actions  taken  by  a  Kalman 
Filter  during  the  crocessing  of  a  measurement  can  be 
obtained  through  the  graphical  interpretation  and  visualiza- 
tion of  these  actions  in  a  two-dimensional  problem. 

From  the  conclusions  of  Section  IV#A,  it  is  advantageous 
to  process  any  group  of  measurements  with  statistically 
independent  noise  components  separately.  In  this  section* 
thuSf  the  processing  of  only  one  measurement  at  a  time  is 
consi  dered. 

1  .   The  Linear  Measurement  Case 

We  will  consiaer  a  dynamic   system   with   two   state 
variables/  and  a  linear  measurement  of  the  form 

z(k)  =  H(k) .x(k)  +  vCk)  = 


=  h^(k).x^(k)  +  h^CkD.x^Ck)  +  v(k) 


with  v(k)  a  zero-meanr  discrete  white  Gaussian  noise  process 
with  E  tv^(k)J  =  r^(k) . 


At  time  t,  a  prediction  is  made  from  previous  esti- 
mates and  knowledge  of  the  plant  and  system  noise  charac- 
teristics. Equations  (4,5)  and  (4.4)  are  used  for   a   linear 
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Figure  8  ;  Comparision  of  processing  policies 
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systernf  (3,7)  and  (3.8)  for  a  general  nonlinear  system. 

The  Kalman  equations  (^.5)/  (4.6)  and  (4.7)  are  now 
applied  to  correct  the  predicted  value  according  to  the  new 
measurement • 

Let's   represent   the    quantities    and    relations 

involved  as  shown  in  Figures  9  and  lOf  where  e   and  e   are  a 

"1      ~2 

pair  of  orthonormal  basis  vectors  and  a   general   vector   is 
represented  by 


^    '-    'l^-l  *  "2^2 


The  following  conventions  and  simplifications  apply 


X*  =  X (k)  = 


X  (k) 


X2Ck) 


X  =  x(k)  = 


=  x(klk-l)  = 


h  = 


"i^    ( k )" 

h 

1 

h^Ck) 

:i 
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tate 


Figure  9  :  Graphical  interpretation  -  Measuremeat 
lines  and  gradient. 


^6 


x;  X* 


Figure  10  ;  Graphical  interpretation  -  Error  ellipse 
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g    =    G(k)    = 


g^Ck) 


Q^^^'^ 


,     ,     _     ,     2    ^       2,1/2 

|g|   =    tg^   +  g2J 


a  =    arctan    (h    /h    ] 


B  =    arctan     ^^^^'S-,^ 


Let's    define    the    function 


<J>  ( X.)     =    h  ,  X,      +    h2  X  2 


(4.26) 


and  call  "measurement  line"  the  locus  of  points   where  <})   is 
constant.   Some  of  these  lines  are  shown  in  Fig  9, 

Using  this  definition  one  sees  that  the  previously 
defined  h  is  the  gradient  of  *  with  resoect  to  x*  or  -^  t 
which  in  this  linear  measurement  case  is  not  a  function  of 
the  states . 

The  prediction  error  covariance  matrix/ 


P(k  !  k-1  )  =  P'  = 


11 


'12 


12 


P22 


is  represented  by  an  error  ellipse  whose  direction  (  e  )   and 
axis   dimensions   (a   -  major  axis?  b  -  minor)  are  given  by^ 
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12] 


tan  2  9  =  2p  '  -^  <  P, '-,   "  °  22  ^   .  .  .  -tt  <  2  9  <  tt 


a  =  (Pj^'^   +  P*22^^^  *  ^12^  ^'"  ^® 


b  =  (D'^^  +  ^22^^^  '  ^W    ®'"  ^ 


(4.27) 
(a. 28) 
(a. 29) 


a.   The  State  Estimate 


Equation  (4,6)  can  be  represented  in  the  form 


X  =  X '  +  g, [resi  dual  1 


(4.30) 


where 


[residual]  =  z(k)  -  H ( k ) . x ( k | k-1 )  = 


=  (j)  (x_*)  -<})(x_')  +  v(k) 


(4.31) 


and  (4.30)  means  that  a  vector  correction  of  magnitude 
jg { . I  res i dual  I  is  made  to  the  predicted  state  in  the  direc- 
tion of  the  vector  a  which  has  an  angle  3  as  shown  in  Fig- 
ure 1 1  . 

Equation  (4,5)  giveSf  as  shown  in  Appendix  A/ 


g  = 


p'  h   +  p'  h 
11  1     12  2 


p'  h"^  +  2p'  h  h   +  p'  h^  +  r 
11  1      12  1  2     22  2 


p'  h   +  p'  h 
12  1      22  2 


'11  \  *    ^^l\S    '    ^22^2  '    '' 


(4.32) 
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(f)(x)  =  <t>{x*) 

(J)(X)  =\(})(X') 


Figure  11  :  Graphical  interpretation  -  Gain  vector, 
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(1)   The  Direction  of  Correction 

The  direction  of  g  is  seen  to  be  indepen- 
dent of  the  measurement  noise^  or,  given  h  and  P(k|k-1)  the 
direction  of  correction  is  already  determined  by 


3  =  arctan  lCPi2^i  +  P22^2^  ^  ^^11^1*  ^12^2^^ 


(a. 33) 


From  Equations  (a, 27)  -  (4.29)  one  obtains 


p: 


12 


11 


P* 


22 


i  (a  -  b).sin  2e 
2 


i  ((a  +  b)  +  (a  -  b)  .cos  2Q I 
2 

i  [(a  +  b)  -  (a  -  b)  .cos  29] 
2 


(a, 3a) 


(4.35) 


(4.36) 


Apolying  these  equations  and  the  definition 
of  a     to  (4.33)  gi  vesr 


tan  3  = 


(a/b  -  l).sin2a  +  ((a/b  +  1)  -  (a/b  -  1)  .cos29]  . t an  a 
(a/b  +  1)  +  (a/b  -  l),cos29  +  (a/b  -1 ) . s i n2 e. t an  a 


(4.37) 


which  can  be  reduced  to 


tan     e   = 


(a/b)  .tane.  [tane.  tana    +13     +     (tana   -    taneJ 
(a/b)  .  [tan  e.  tana    +    11     -    tane.  (tana    -    tanSJ 


(4.3a) 
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This  equation  shows  that  the  direction  of 
the  correction  generated  by  the  Kalman  Filter/  in  this  two- 
state  oroblem^  is  a  function  only  of  the  direction  (  a  )  of 
the  gradient  (h^)  and  of  the  alignment  (0)  and  the  ratio 
(a/b)  of  the  prediction  error  ellipse. 

Three  special  cases  of  interest  occur  when 
(i)  a/b  =  1/  (ii)  b  =  0,  and  (iii)  e=a  • 

Case  ( i )  -  If  a/b  =  t  the  error  "ellipse" 
is  actually  a  circle  meaning  that  the  prediction  has  the 
same  uncertainty  in  any  direction.  No  preferable  direction 
pre-exists  and  the  correction  is  in  the  direction  of  the 
gradient  of  the  measurement  function^  as  can  be  seen  from 
(a. 37) 

tan  B  =  (0  +  2,tana)  /  (2  f  0)  =  tana 


or 


=  a  -f  n  IT 


Case  ( i  i )  -  If  b  =  0  the  error  ellipse  is  a 
line  meaning  that  the  prediction  is  exact  along  the  minor 
axis  and  that  only  corrections  along  the  major  axis  are 
necessary.  As  can  be  anticipated/  the  direction  of  correc- 
tion will  be  the  direction  of  the  major  axis#  as  can  be  seen 
by  taking  the  limit  as  b->-0  of  Eg,  (^.30),  which  yields 

tan  3  =  tan  9 
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or 

3  =  e  +  n  IT 

For  any  other  intermediate  value  of  a/b  between  1  and   »   , 
the  value  of  3   will  be  between  a   and   9   ( ■♦'  n  tt)  . 

Case  ( i  i  i  )  -  When  9  =  a  the  gradient  is 
colinear  with  the  major  axis  of  the  prediction  error  ellipse 
an6f     from  (a. 38), 

tan3  =  [(a/b) .tan9. (tan^e +  1)  +  01  /  [ia/b).itar?e^    1)  -  0] 
or 


tan3   =  tan 


and 


=  9  +  nTr=  a+n-iT 


Thus  the  direction  of  the  correction  is 
along  the  major  axis  of  the  ellipse  which  is  colinear  with 
the  gradient  of  the  measurement  function. 

(2)   The  Amount  of  Correction 

The     magnitude     of  the     correction 

(  {  resi  dua  1  I  .  I  g  l| )  can  be  better  seen  if  one  decomposes  g  into 

components  m  and  n,  along  the  gradient  and  the   tangent   to 

the  measurement  function,  as  shown  in  Fig,  11, 


m  =  g,cosa+  g  .sina 


and  from  the  definition  of  a. 


(a. 39) 
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COS  a  =    h  .  /  |]i  I     t 


sina=    h/jh_j 


(a. 40) 


Using    Equations     (a. 32),     (a. 39)    and    (4.aO), 


m    = 


p'     h^     t    p'     h    h       +    p'     h    h       +    p'     h 


11     1 


12     1    2 


12     1    2 


p'     h^     f    2p'     h    h       +    p'     h^     +    r 
11   1  12      1      2  22      2 


22      2 
2 


ly 


1     +Y 


ly 


ca.ai) 


where 


2 
Y     - 


/    (p'    h       +    ao'     h    h       +    p'    h    ) 
11     1  12      1      2  22      2 


(4.42) 


Consider  some  special  cases.  How  would   the 

new   estimate   be   computed  if  a  noise-free  measurement  were 

obtained  and  the  Kalman  filter  were   aware   of   this   fact» 
i.e.,  r  =  0? 

Note  from  Equations  (4.41)  and  (4.42)   that 
when  the  measurement  is  noise-free,  r  =  0,  y  =  0  and 


=  1  /  jhf  , 


a  constant. 


No  matter  what  P(k!k-1)  is,  the  projection 
of  the  correction  into  the  gradient  is  a  constant  if  r  =  0, 
or,  the  tip  of  the  correction  vector  follows  a  line  normal 
to  the  gradient  when  P(k|k-1)  is  varied. 

Figure  12  shows  the  correction  made  to  the 
predicted   state   when  b  =  0.  In  Appendix  C  it  is  shown  that 
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for  the  present  case  of  a  linear  measurement  function  the 
vector  OA,  in  the  direction  of  h  and  length  ( res i dua 11 /  |  h J 
has  its  end  point  A  on  the  line 

(J)  (x_)  =  «|'(x_')  f  residual  =  (f)  ( x.*  )  +  (j)  (^* )  ♦  v  -  <|)  ( x^' ) 


or 


<i>  (x)  =  <}>(x*)  +  V  =  2 


whichf  for  r  =  0  (i.e.  v  =  0)  passes  through  the  point  x_      = 
x_*. 

Since  the  correction  has  the  direction  of 
the  prediction  error  ellipse  (line)  and  has  OA  as  its  pro- 
jection into  the  vector  h_»  it  can  be  seen  from  Figure  12 
that  in  this  case  of  b  =  0  and  r  =  0  the  estimate  will  coin- 
cide with  the  true  state. 

Figure  13  shows  the  correction  for  a  gen- 
eral error  ellipse  but  still  with  r  =  0.  It  is  seen  that  the 
projection  into  h  is  the  same  as  before.  The  direction  of 
correction  3  is  somewhere  between  a  and  6  • 

A  general  correction  for  a  noisy  measure- 
ment is  shown  in  Figure  l^r  with  different  scaling  than  the 
previous  figures.  The  direction  of  correction  3  is  first 
determined  and  suppose  it  is  as  shown.  The  projection  of  the 
correction  into  h^  would  be  OA  if  r  were  zero?  for  r  ^  0  the 
reduction  factor  is  applied  and  suppose  the  projection  is 
vector  OB  as  shown  in  the  figure.  Taking  a  normal  to  h   from 
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(j)(x)  =  z 


error  ellipse  (line) 


perfect  estimation 


^.  res  idual 


length  =  residual/  [h ( 


Figure  12  :  Linear  measurement  -  Perfect  estimation 
for  b=0  and  r=0 . 
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<t)  (x)  =  z 
(t)(x)  =\(f)(x') 


^  ^.residual 


residual  /  |h| 


Figure  13  :  Linear  measurement  -  r=0 
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point  Bf  the  point  C  where  it  cuts  the  direction  of  correc- 
tion determines  the  final  correction  OC  and  the  new  esti- 
mate • 

b.   The  Estimation  Error  Covariance  Matrix 

The  effects  on  the  estimation  error  covariance 
matrix  can  be  seen  in  the  following  way.  From  Equation 
(a. 7), 

P(k)  =  II  -  G(k)H(k)lP(k|k-l) 

or 


11 


12 


12 


22 


1  -  g  h 
1  1 


g  h 
2  1 


-  q  h 
1  2 


1  -  g  h 
2  2 


P'        P 
11       12 


12 


P' 

22 


Using  the  values  of  g,  and  q^      from   (^,32)   one 


finds    that 


Pll=     ^^2^Pl1  ^22     "    Pl2  )     +    Pil    --^l     /    c 


(a.asa) 


^12=     '\'2'-il    -    ^U^22^     '    ^12^'^     '    ' 


p      =     Ih^Cp'    p' 
22 


-d'^)    +p'     r2]     /c 
1      11      22  12  22 


(a.a3b) 
(a.a3c) 


where 


C       =     [p',   h^    +    2D'hh       fp'h^    +    r^] 
11     1  12      1     2  22      2 


(a.a3d) 
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<i>{x)   =  (t)(x*) 

^  (X*)   +  V  =  2 


residual  /  |_h| 


i 


Figure  14  :  Linear  measurement  -  General  correction, 
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And  the  estimation  error  ellipse  has  the   orien- 


tat  i  on 


tan  26  =  ^Pi2'^  ^^n  "  ^22  ^  ' 


(h: 


^^^2^^'l2  ■  ^U^22^  '    '^U^ 


-  p'  )r 
22 


(^.a^) 


The  limiting  situations  occur  in  the  special 
cases  when  (i)  r  ->  oo  and  (ii)  r  =  0, 

Case  ( i )  -  When  r  is  relatively  large*  little 
information  is  provided  to  the  filter;  the  estimate  is 
essentially  the  same  as  the  prediction*  and  the  estimation 
error  ellipse  is  essentially  the  same  as  the  prediction 
error  ellipse.  This  can  be  seen  from  Eguations  i^.'iia) 
(a.aSd)  f    as  r  ->  00  . 

Case  (ii)  -  when  r  =  0*  all  errors  in  the  pred- 
iction* in  the  direction  of  the  gradient  of  the  measurement 
function*  are  corrected  and  the  estimation  error  ellipse 
becomes  a  line  colinear  with  the  measurement  line.  As  shown 
below*  this  result  is  obtained  from  Eguations  (4,43)  and 
(4.44)  when  r  is  made  egual  zero. 


From  (4.a3a)  -  (4.43d)*  p   ,p    is  easily   shown 

11   22 
2 
to   be  egual  to  p    when  r  =  0*  indicating  that  P(k)  becomes 

12 

singular.  From  (4.44), 


tan  26  =  2h^h  ^  /     ( h^  "  ^  2^  "^ 
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=  2(h  /h  )  /  (1  -  (h  /h  )2  ]  = 

2   1  2   1 

2 

=  2  tan  a  /  (1  -  tan  a)  =  tan  2a 

From  analysis  of  the  individual  terms  of  (^.43a) 

(a,a3d)»  this  result  is  to  be  interpreted  as 

2  6   =  2a   +  TT 


OP 


6   =  a   f   tt/2 


When  the  minor  axis  of  the  prediction  ellipse  is 

2 
already   zeror   as  in  Figure  12f  ^,e,^    d'   =  p*  p'  f  and  the 

measurement  is  noiseless^  then  P(k)  =  0   and   complete   cer- 
tainty about  the  system  state  is  obtained. 


These  two  last  oaraqraphs  show  a  very  interest- 
ing situation:  if  two  perfect  measurements  are  made  on  a 
two-state  plant  and  the  filter  is  aware  of  this  (i.e.  r  =  0 
is  used  in  the  Kalman  equations)*  then  an  exact  estimate  is 
obtained  ---  the  estimation  error  covariance  matrix  col- 
lapses into  the  zero  matrix. 

For  a  general  situation  with  normal  measurement 
noise  values*  the  estimation  error  ellipse  is  rotated  toward 
the  measurement  line*  away  from  the  gradient*  with  5  occu- 
pying a  position  between  e  and  a  +  Tr/2. 
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The  general  expressions  for  the  lengths  of  the 
major  and  the  minor  axis  of  the  estimation  error  ellipse  can 
be  obtained  by   the   application   of   Equations   (^.28)   and 

(a. 29)  to  (a.asa)  -  (a,a3c). 

c.   The  Important  Points 

Summarizing  the  principal  points  seen  up  to  now* 

---  the  direction  of  correction  in  the  estimate  is 
defined  at  the  end  of  the  prediction  phase  and  doesn't 
depend  on  the  measurements. 

---  if  the  measurement  is  noiseless  and  r  =  0  is  used 
in  the  Kalman  filter  equations^  the  projection  of  the 
correction  on  the  gradient  of  the  measurement  function  is 
(residual  /  (jlp'  independent  of  P(k{k-1).  The  estimation 
error  covariance  matrix  becomes  singular  indicatingr  in  a 
two-dimensional  problem*  that  the  estimation  error  ellipse 
becomes  a  line  ---  a  point  when  the  prediction  error  ellipse 
is  al ready  aline. 

---  if  the  measurement  is  noisy*  the  projection  of   the 

correction   on   the   gradient   is  multiplied  by  the  reducing 

2 
factor  (1  /  (1  +  Y  )],   The  error  ellipse  is  rotated   toward 

the  measurement  line*  away  from  the  gradient  of  the  measure- 
ment f unct  i  on  . 

---  if  the  measurement  is  too  noisy*  the  reducing  fac- 
tor  tends   to   zero   and   the  estimation  is  the  same  as  the 
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predi  ct  i  on. 

2.   The  Nonlinear  Measurement  Case 

Assume  the  same  two-dimensional  system  as  before  but 
with  nonlinear  observations  of  the  form 


2(k)  =  h(x(k),t  )  +  v(k) 

—        K 


a.   The  Extended  Kalman  Filter  Approach 

The  prediction  phase  is  the  same  as  before  but 
the  estimation  equations  are  now  (3,18),  (3,19)  and  (3,21) 
for  an  Extended  Kalman  Filter, 

Let's,  again  define  a  measurement  function  by 

^(x^)  =  h(x^) 

and  redefine  Ji  as  the  gradient  of  9   with  respect  to  x^. 

Now  the  components  of  h  are  no  longer  constants 
but  depend  on  x_.  At  the  predicted  point  these  components 
will  be 


h   =4-  h(i(k|k-l))    , 


3X 


h    =-^  h(i(k;k-i)) 

3x«    — 


Equations  (a. 27)  through  (n,HH)       are   all   still 


valid. 


Let's  consider  two  typical  measurement  functions 
and  observe  graohically  how  the  new  estimate  is  obtained. 
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Case  (i)  -  assume  the  measurement  function  is  of 


the  form 


((>(_x)  =  arc  tan  Kx   -  a  )  /  (x   -  ^9  ^  ^ 


(a. 45) 


The  measurement  lines  will  be  as  shown  in  Figure  15, 

Suppose  now  a  noiseless  measurement  is  obtained 
and  is  to  be  processed.  Fig,  16  shows  the  situation  when  the 
prediction  error  ellipse  is  a  line  (b  =  0),  In  Appendix  C  it 
is  shown  that. for  this  special  measurement  function  the  vec- 
tor OA  of  length  [resi dual  J / 1 h I  will  end  somewhere  before  a 
point  B  on  the  curve  (\)(.x_)  ~  z.  Supposing  its  length  is  as 
shown/  it  can  be  seen  that  an  imperfect  estimate  is 
obtained.  Nevertheless  the  estimation  error  covariance 
matrix*  given  also  by  Equations  (4.43)f  will  collapse  into 
the  zero  matrix  indicatingr  falsely*  that  a  perfect  estimate 
was  obtained. 

Figure  17  shows  a  typical  situation  for  a  gen- 
eral prediction  error  and  a  noise-free  measurement.  The  tip 
of  the  correction  vector  can  be  seen  to  follow  the  line  A 
A'  parallel  to  the  line  (j)(x^)  =  c()  ( x^' )  instead  of  being  along 
the  line  ^(.x_)  =  (jjCx^*).  As  seen  in  the  previous  section, 
the  estimation  error  ellipse  will  erroneously  be  determined 
as  a  line  also  along  A  -  A',  indicating  falsely  the  direc- 
tion of  the  true  state. 


Case  (ii)  -  assume  the  measurement   function   is 


of  the  form 


loa 


^2 


e 
-2 


C{)  (X)  =  (J,  (x*) 


=  (|)(x') 


Figure  15  :  Nonlinear  measurement  lines  I 
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^.  residual 


. -  residual  /Ihl 


(J)  (x)  =  cf)  (x*)  =  z 


J(x)  =  (t)(x') 


Figure  16  :  Nonlinear  measurement  I  -  Imperfect 
correction  for  b=0  and  r=0 . 
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^.  residual 


residual  /  |h| 


(|)(x)  =  (1)(X*)  + 


V  =  z 


({)(X)  =  (t)(x') 


Figure  17  :  Nonlinear  measurement  I 
correction  for  r=0 . 


-  General 
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^( 


X )  =  ( X   -  a 


2  2 


1/2 


(a.ab) 


The  measurement  lines  are  as  shown  in  Figure  18.  This  same 
figure  shows  the  new  estimate  obtained  when  the  measurement 
is  noiseless  (r  =  0)  and  the  prediction  error  ellipse  is  a 
line  (b  =  0),  The  same  observations  of  Case(i)  are  valid. 
As  shown  in  Appendix  C  the  point  A  is  on  the  curve  (j)(x)  =  z 
for  this  special  measurement  function. 

These  figures  have  their  values  mostly  because 
they  showf  in  a  simple  wayr  the  large  errors  that  can  be 
made  in  the  processing  of  nonlinear  measurements  by  an 
Extended  Kalman  Filter.  Two  types  of  errors  are  generated: 
one  in  the  determination  of  the  new  state  estimate;  the 
other  is  seen  in  the  fact  that  the  estimation  error  covari- 
ance  matrix  P  becomes  a  poor  approximation  to  the  true  error 
covariance^  making  the  filter  non-optimal.  These  difficul- 
ties are  a  function  of  the  non 1 i near i t i es  involved  andf 
apparentlyr  are  more  significant  when  the  predicted  values 
are  poor  and  when  the  measurement  noise  is  assumed  to  be 
small  when  evaluating  the  extended  Kalman  filter  equations. 

b.   The  Iterative  Approach 

A  way  of  correcting*  in  part*  these  deficiencies 
is  the  adoption  of  iterative  procedures.  The  first  improve- 
ments can  be  obtained  if*  after  the  Extended  Kalman  Filter 
has  been  applied  to  correct  the  prediction*  the  filtering 
process  is  repeated  but  with  the  linearizations   made   about 
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(|)(x)  =  <t>{x*)    =   z 


Figure  18  :  Nonlinear  measurement  II  -  Imperfect 
correction  for  b=0  and  r=0 . 
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the  just  obtained  estimate^  hopefully  a  better  approximation 
to  the  true  state  than  the  crude  prediction. 

Analytically  the  process  consists  of  replacing 
the  estimation  equations  (3.18)*  (3.19)  and  (3.21)  by  the 
iterative  equations 


X.    =  j((k|k-l)  t  G  [z(k)  -  _h(x^(k|k-l))l 
i+1  i 


(a. a?) 


G.  =  P(k|k-l)Hj[H.  P(k|k-1)H!'  f  R]  ^ 


(a.as) 


H  .  =  -2-  h(x.  ) 


3x 


—1 


f  i=0,  1»2/...,  f 


where 


X 

— o 


=  x(k!k-l)    , 


X   =  x(k) 

-f    - 


P(k)  =  tl  -  G^H  JP(k!k-l) 

f  f 


(a.a9) 


For  our  two-dimensional  problem  Equation   (^.47) 
can  be  expressed  as 


=   X '  +  q  ,  [resi  dual ] 
i+1    -     ^i 


(a. 50) 


where 


residual  =  z  -  h(x')  =   <J)  ( x  * )  -  (}>  ( x  '  )  +  v(k)     (4.51) 


Adding  and  subtracting   ((>  ( x_  )  to  the  residuals 

i 
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X    =  x'  +  g  .  fz  -  rf,(x  )  •♦•  (J)(x  )  -  (}>(x')] 
-i+1   -      i  *-i       -i     '^- 


or 


x.^  =  x'  +  g_..[z    -  <j)(x.)l  -  q..[(J.(x')  -(l)(x.)]   (4.52) 


These  equations  can  be  applied  until  there  is  no 
significant  difference  between  consecutive  estimates^  and  2 
or  6    iterations  are  normally  enough. 

The  effect  of   these   equations   can   be   easily 

visualized   for   the   simple   case   of   b  =  0  and  r  =  0.  The 

sequence  of  Figures  19  -  21  show  this  effect  for  one  of   the 

special  measurement  functions  already  studied.   In  Figure  19 

the  first  estimate  is  obtained  as  done  previously  in   Figure 

16.   The   vector   OA,   which   is  smaller  than  OB  as  shown  in 

Appendix  Cr  represents  the  projection  of  the   correction   on 

the  gradient  h  ,   A  normal  to  h   from  point  A  cuts  the  pred- 
—o  — o 

iction  error  line  (b  =  0)  at  the  first  estimate  x  . 

In  Figure  20  the  terms  of  Equation  (4,52)  are 
shown  for  i  =  1.  Two  correction  terms  are  applied  to  the 
prediction  x.' ,  both  in  the  direction  of  g^  but  with  dif- 
ferent residuals.  The  direction  of  g  is  again  along  the 
prediction  error  line  (b  =  0). 

The  gradient  of  the  measurement  function  at  the 
first  estimate  is  Ji  .  The  projection  of  the  first  correc- 
tion (g.,lz  -  (i>(x..)l)  into  h^^  is  smaller  than  CEr  say  CO. 
Taking   the   normal   at  D  one  gets  the  first  correction  term 


HI 


irs  t 
estimate 


actual 
state 


(})  (x)  =  (j)  (x*)  =  z 


(x)  =  (})(x')  1     11 


/ 


Figure  19  :  Iterative  process  -  First  estimate. 
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CH. 


The  projection  of  the   second   correction   term, 

(-g^  I  <t>  (.)t_*  )       -  <t>  (i_    )})f  into  h   is  smaller  than  CG/  say  CF, 
i  i  —1 

The  normal  from  point  F   determines   the   second   correction 
term  (IC). 


Adding  the  three  vectors  x^*  f  CH  and  IC  one  has 
the  second  estimate,  shown  in  Figure  21/  and  the  next  itera- 
tion can  be  processed. 


If  in  Equation  (^.^7)  one  modifies  the   residual 


by 


2(k)  -  h(x(klk-n)  =  z(k)  -  fhCx  )  + 


+  H  .(x(k|k-n 
i   — 


X  )  +  .  .  .] 
—  i 


=  z(k)  -  h(S  )  -  H  ,  [5(k|k-l)  -  i  ] 

—       —    —L  i   —  — i 


then  Equation  (^.^7)  becomes 


i  ,      =    5(k|k-l)    +    G     Izik)    -    h(x    )    -    H     [xCklk-n    -    X    ]] 
—  1+1      -  1—  i  X  —  — i 

(a. 53) 


This  variation  in  the  residual  is  very  small  to 
affect  considerably  the  estimates  but  in  till  it  is  shown 
that  its  use  produces  a  better  agreement  between  the  calcu- 
lated estimation  error  covariance  matrix  and  the  true  error 
covariance,  when  processing  a  noisy  measurement. 
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^2 


y 


|(|)(x*)-(t)(x^)|/|hj 


(|)(X)   =  (()(X^) 


(J)(X)  =  (|)(X*)  =  z 


^1 


(i)(x)  =  (f)(x') 


Figure  20  :  Iterative  process  -  Correction  terms. 
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second 
estimate 

(J)  (x)  =  (|)  (x^) 


=  (f)(x*) 


(j)  (x)  =    (t>  (X2) 


(J)(X)  =  (J)(X') 


Figure  21  :  Iterative  process  -  Second  estimate 
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3.   Summary 

The  direct  application  of  these  graphical  interpre- 
tations to  a  general  problem  is  impractical.  Nevertheless 
the  basic  points  observed  here  can  be  generalized  for  Kalman 
Fi 1 ters: 

---  the  ••direction"  of  correction  in  the  estimation 
is  defined  by  the  predicted  quantities  and  doesn't  depend  on 
the  measurements. 

---  the  amount  of  correction  in  an  estimation,  along 
the  direction  defined,  is  the  result  of  a  weighting  process 
between  the  uncertainty  in  the  measurement  and  the  confi- 
dence in  the  predicted  values.  This  correction  is  a  function 

2 
of  the  factor  [1  /  (1  +y  )1  where,   in   general,   for   a   n- 

dimensional  system 

n   n 

Y^  =   r^  /.^,  \        d'  h.h  (a. 51) 

1-1  j=l    ij  1  J 

---  the  errors  in  the  estimate  from  the  processing 
of  a  nonlinear  measurement  can  be  large  and,  most  important, 
the  estimation  error  covariance  matrix  may  become  a  very 
poor  approximation  of  the  true  error  covariance,  forcing  the 
filter  to  inaccurately  weight  future  measurements.  This 
problem  depends  on  the  nonl i near i t i es  involved  and  seems  to 
be  worse  when  the  predicted  values  are  poor  and/or  the  meas- 
urements are  considered  by  the  filter  to  have  relatively  low 
noi  se. 
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---  the  errors  indicated  above  may  be  reduced  with 
the  use  of  the  iterative  equations  (a.aS),  (a.a9)  and  (a.a7) 
or  ('4^.53)r  instead  of  the  standard  approach  given  by  the 
Equations  ( 3. 18) r (3. 19)  and  (3.21)^  at  the  expense  of 
greater  computing  effort. 
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V.   COMPUTER  SIMULATION 


A.   GENERAL  CONSIDERATIONS 

A  computer  simulation  was  implemented  in  a  Unix  Operat- 
ing System  running  on  a  POP  11/50  with  the  FP  IIB  floating 
point  processor^  located  at  room  506  of  Spanagel  Hal  1 r  NPS. 
The  whole  program  is  interactive  and  is  expected  to  be 
self-explanatory  to  a  user.  Its  block  diagrams  and  coding 
are    presented  in  Appendix  E. 

The  graphical  outputs  were  generated  on  a  Tektronix 
aOl^-l  terminal.  The  presentation  of  the  results  of  the 
various  simulations  is  simolified  by  the  adoption  of  the 
following  conventions. 

The  metric  system  was  used  in  all  calculations  and  out- 
puts. The  axes  are  dimensioned  in  kilometers.  Heading  and 
bearing  angles^  expressed  in  radiansr  are  measured  counter- 
clockwise from  the  positive  X-axis. 

Figure  22  shows  a  typical  plot.  Points  of  the  true  track 
are  represented  by  the  sysmbol  "+"  and  interconnected  by  a 
solid  line.  Buoys  are  represented  by  a  distinctive  boxed 
symbol.  The  portion  of  the  true  track  that  is  in  the  range 
of  a  buoy  is  bounded^  if  necessary^  by  dash-aotted  lines. 
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The  signals  received  by  the  sonobuoys  are  assumed  to  be 
sent  to  a  ship  or  aircraft  where  they  are  processed  by  spe- 
cialized equipment  to  generate  certain  measurements.  The 
major  characteristics  of  these  measurements  are  printed  on 
top  of  each  figure  and  obey  a  code  whose  general  form  is 

L   Nl   N2   N3   S 

where 

L   -  a  capital  letter  indicating  a  type  of  measurement. 

B   stands   for   bearing?   F   for   frequency?  D  for  frequency 

difference;  R  for  frequency  ratio;  and  T  for  time  delay. 

Nl  -  an  optional  number  indicatingf  when  necessaryr  the 
buoy  that  generates  the  measurement.  If  this  number  has  two 
digits  it  identifies  two  buoysr  thus^  T  12  indicates  a  time 
delay  measurement  between  buoys  1  and  2, 

N2  -  a  number  indicating  how  many  seconds  after  the 
start  of  the  problem  this  measurement  will  be  first 
obtained.  It  is  also  optional. 

N3  -  a  number  indicating  the  period  between  consecutive 
measurements. 

S  -  a  number  indicating  the  standard  deviation  of  the 
simulated  zero-mean  Gaussian  noise  that  is  associated  with 
this  measurement/  followed  by  an  ootional  small  letter  indi- 
cating dimension:  d  for  degrees^  h  for  Hertz  and  s  for 
seconds • 
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The  interpretation  of  the  code  presented  in  Figure  ?2  is 
then  that  both  buoys  1  and  2  provide  bearing  and  frequency 
measurements  with  a  period  of  100  seconds  between  consecu- 
tive set  of  measurements.  Buoy  1  provides  its  first  set  of 
measurements  100  seconds  after  the  start  of  the  problem; 
buoy  2r  perhaps  because  of  its  limited  ranger  provides  its 
first  set  of  measurements  only  after  500  seconds  of  the 
start  of  the  problem.  The  standard  deviation  of  the  measure- 
ment noise  is  5  degrees  for  the  bearings  and  0.04  Hertz  for 
the  frequencies. 

The  filter  receives  the  measurements  for  processing.  It 
is  given  some  initial  conditions  and  the  initial  position  is 
marked  in  the  figures  by  the  first  diamond  with  the  letters 
I. P.  on  top. 

As  the  target  moves  and  the  measurements  are  simulatedr 
the  filter  generates  estimates  of  the  target  parameters 
which  are  represented  by  small  diamonds^  interconnected  by 
dashed  lines  to  represent  the  estimated  path.  The  diamonds 
and  the  "+"  symbols  correspond  to  the  same  time  points. 

Associated  with  each  estimated  point  the  filter  provides 
an  estimation  error  covariance  matrix  which  indicates  the 
amount  of  confidence  or^  conversely^  uncertainty  that  should 
be  given  to  this  estimate.  An  interesting  way  to  show  this 
in  a  trajectory  plot  like  Figure  22  is  by  plotting  the  one- 
sigma  error  ellipses  associated  with  each  estimated  point. 
These  error  ellipses  have  their  orientation  and  axis   dimen- 
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sions  given  by  Equations  (^,27)  -  (a, 29),  and  can  be  seen  in 
Figure  23. 

These  figures  were  obtained  by  running  the  simulation 
once.  They  may  represent  a  typical  real-time  result. 
Nevertheless  they  do  not  give  the  expected  or  average 
behavior  of  the  filter  to  the  situation  simulated.  To  obtain 
this  average  behavior  many  Monte  Carlo  runs  have  to  be  exe- 
cuted and  their  individual  results  processed  to  generate  the 
sample  statistics.  The  sample  statistics  formulas  for 
obtaining  the  sample  means  and  covariances  can  be  found,  for 
example,  in  page  86  of  Ref.  [2J  • 

In  this  work  all  the  statistical  results  were  obtained 
from  200  MonteCarlo  runs.  Figure  2^  shows  the  average 
behavior  of  the  filter  after  these  runs.  Now  the  small  dia- 
monds represent  the  sample  mean  estimated  position  and  the 
ellipses  show  the  spreading  of  the  results  about  the  mean 
values.  Note  that  at  the  initial  position  point,  since  in 
all  runs  it  was  the  same,  the  ellipse  is  the  point  itself. 

To  observe  the  errors  in  the  estimation  of  the  speed  and 
heading  of  the  target  different  plots  have  to  be  used.  Fig- 
ure 25  shows  the  sample  mean  error  in  the  estimation  of 
speed  as  a  function  of  time.  Figure  26  shows  the  sample 
mean  error  in  the  estimation  of  the  heading  of  the  target. 

Other  types  of  plots  or  outputs  could  be  generated  from 
the  results  given  by  the  program  simulation,  like  the  sample 
covariance  between  estimated  velocity   and   bearing,   if   it 
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ever  becomes  necessary. 


B.   SELtCTED  RESULTS 

As  has  been  seen  in  this  workr  the  number  of  variables 
in  the  passive  tracking  problem  is  very  large.  Some  specify 
the  true  target  track  and  the  true  position  and  characteris- 
tics of  the  sonobuoys;  others  define  the  measurement 
scheduling  and  noise  chara<:t  er  i  st  i  cs  .  The  nonlinear  filter 
that  processes  the  measurements  has  its  own  variables  and 
must  also  assume  values  for  the  other  variables^  which  may 
be  different  from  the  true  ones. 

For  any  combination  of  these  variables  a  new  problem  is 
formed  and  its  simulation  may  generate  a  number  of  interest- 
ing conclusions.  The  results  presented  are  representative  of 
those  obtained  from  the  execution  of  the  simulation  program 
described  in  Appendix  E.  In  most  of  the  results  presented  a 
certain  reasonable  initial  condition  for  the  filter  was 
assumed;  the  position  of  the  buoys  was  assumed  constant  and 
known;  the  ranges  within  which  reliable  bearing,  frequency 
and  time  delay  measurements  could  be  obtained  were  assumed 
e(5ual  . 

1.   Non-maneuvering  Target 

a.   One  Buoy  in  Action  at  a  Time 

Figure  27  shows  the  behavior  of  the  filter  in   a 
situation   where  only  one  buoy  is  in  contact  with  the  target 
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at  any  instant.  The  range  and  position  of  the  buoys  were 
pre-arranged  to  satisfy  this  requirement.  Target  speed  is 
5.0  m/s?  its  heading  is  315  degrees  measured  from  the  posi- 
tive X-axis.  The  initial  assumptions  of  the  filter  include 
6.0  m/s  for  speed  and  325  degrees  for  heading.  The  measure- 
ments were  processed  as  suggested  in  [21.  Frequency  and 
bearing  measurements  are  available  every  100  seconds  and 
processed  simultaneously?  standard  deviations  of  5  degrees 
and  0.04  Hertz  were  assumed  for  the  noise  components. 

The  frequency  measurement  provides^  indirectly^ 
some  range  information  and*  this  way/  complements  the  bear- 
ing measurement.  This  range  information  can  also  be  obtained 
by  judiciously  positioning  the  next  buoys  or  by  having  two 
or  more  buoys  in  contact  with  the  target  at  one  time. 

In  Figure  28  the  frequency  measurement  was  elim- 
inated and  the  simple  bearing  measurement  allowed  the  esti- 
mator to  track  the  non-maneuvering  target  almost  as  well  as 
before^  because  of  the  special  position  of  the  buoys.  Note 
that  the  ellipses  tend  to  align  with  the  measurement  lines. 
This  can  be  explained  from  the  discussions  of  Section  IVrB 
and  Appendi  x  D. 

There  are  many  ways  to  improve  this  tracking.  A 
very  small  improvement  can  be  obtained  if  less  noisy  or  more 
frequent  frequency  measurements  could  be  obtained.  To  get 
more  precise  frequency  indication  one  needs  better  resolu- 
tion from  the  FFT  processors;  to   obtain   better   resolution 
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one  needs  longer  records  whichr  in  turn*  introduce  more 
errors  due  to  target-sensor  geometry  variations  during 
record  time.  To  reduce  the  period  between  consecutive 
independent  freguency  indications  one  has  to  sacrifice  reso- 
lution but  this  would  reduce  even  more  the  effect  of  the 
freguency  measurement  on  tracking  accuracy.  Besides*  an 
emitted  freguency  of  300  Hertz  from  a  target  with  a  relative 
velocity  of  5  knots  toward  a  buoy  is  shifted  by  only  about 
0.5  Hertz  and  it  can  be  seen  that  not  much  resolution  can  be 
spared. 

The  bearing  information  is  basically  continuous 
and  it  seems  reasonable  to  admit  that  only  economical  limi- 
tations exist  for  obtaining  less  noisy  or,  egui val ent 1 y f 
more  freguent  bearing  measurements.  Figure  29  shows  the 
effect  of  processing  only  bearing  measurements*  but  more 
frequently  (at  25-second  intervals).  In  Figure  30  the  stan- 
dard deviation  of  the  measurement  noise  was  reduced  to  1 
degree.  Note  that  the  alignment  of  the  ellipses  with  the 
measurement  lines  is  more  pronounced  with  less  noisy  meas- 
urements. Mote  also  that  the  error  in  range  while  the  first 
buoy  is  in  contact  with  the  target  is  not  corrected  by 
improving  the  accuracy  of  the  bearing  measurements. 

Figures  31  and  32  show  the  mean  errors  in 
estimating  the  velocity  and  bearing  of  the  target.  The  sym- 
bol "•♦■"  is  associated  with  Fig  27#  the  triangles  with  Fig  28 
and  the  diamonds  with  Fig  29. 
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Since  the  contribution  of  the  frequency  measure- 
ments is  so  small  compared  to  that  of  the  bearing  indica- 
tions in  this  casef  the  partitioning  of  the  measurements  has 
very  little  effect  on  tracking  accuracy.  Also  the  iterative 
techniques  suggested  in  Section  IV/B,2  have  no  marked  effect 
on  the  processing  of  the  frequency  measurements  for  the  same 
reason. 

b«   Two  Buoys  in  Action  at  a  Time 

The  placement  of  another  buoy  in  action  can 
greatly  improve  the  tracking  accuracy^  as  shown  in  the  next 
figures.  The  buoys  are  now  arranged  so  that  two  buoys* 
instead  of  one*  are  in  contact  with  the  target  at  any  time. 
The  scaling  is  changed  so  that  the  estimation  errors  may  be 
better  appreciated. 

In  Figure  33  the  added  buoy  is  a  Lofar  buoy  that 
can  only  provide  frequency  measurements.  The  two  frequency 
and  one  bearing  measurements  available  every  100  seconds 
were  processed  simultaneously.  Note  the  improvement  in  accu- 
racy when  buoys  3  and  4  gain  contact  with  the  target  in 
replacement  of  buoys  1  and  2.  This  also  is  explained  from 
the  discussions  of  Section  IV^B. 

In  Figure  3^  all  buoys  contribute  with  frequency 
and  bearing  measurements  which  are  still  processed  simul- 
taneously every  100  seconds.  The  position  of  the  buoys*  in  a 
line   perpendicular   to   the   true  track*  is  not  a  very  good 
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position  for  two  Difar  buoys.    Nevertheless   the   resulting 
tracking  accuracy  is  good  for  this  non-maneuvering  target. 

In  Figure  35  the  two  frequency  measurement s^ 
which  require  a  five-state  plant  for  their  processing^  were 
replaced  by  a  frequency  difference  measurement  with  O.OU 
Hertz  of  standard  deviation  in  its  noise  component»  which 
requires  only  a  four-state  plant*  as  discussed  in  Section 
II/D#2.  The  results  are  basically  the  same*  mainly  because 
they  are  mostly  dependent  on  the  bearing  measurement  and  not 
on  the  frequency  information.  The  use  of  the  frequency  ratio 
measurement  also  gave  the  same  kind  of  results  when  a  very 
low  noise  measurement*  with  about  0.0001  units  of  standard 
deviation*  was  considered. 

In  Figure  36  the  frequency  measurements  were 
eliminated  and  only  bearing  measurements  were  processed. 
Because  of  the  bad  position  of  the  buoys*  on  a  line  perpen- 
dicular to  the  true  track*  the  accuracy  of  the  tracking 
along  that  line  is  not  good.  The  ellipses  show  a  reasonable 
spreading  of  the  estimates  along  perpendiculars  to  the  true 
track.  Figure  37  shows  the  same  plot*  with  different  scal- 
ing* for  the  first  1000  seconds  of  the  path. 

For  different  positions  of  the  buoys*  especially 
if  one  moves  the  buoys  off  the  perpendicular  to  the  true 
track*  better  tracking  can  be  obtained  as  shown  in  Figure 
38.  The  tendency  of  the  ellipses  to  align  with  the  measure- 
ment lines  explains  the  improvement  obtained. 
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By  reducing  the  period  between  consecutive  bear- 
ing measurements  and  processing  them  separately  one  may 
obtain  the  accuracy  shown  in  Figure  39.  The  ellipses  are 
like  the  ones  of  Figure  58  but  with  reduced  dimensions. 

Good  results  can  also  be  obtained  if  one 
includes  the  time  delay  measurements.  Figure  ^0  shows  the 
result  when  they  are  applied  to  this  situation.  Notice  that 
the  ellipses  are  very  small  reflecting  the  fact  that  very 
good  information  about  the  position  ^of  the  target  is  given 
by  a  low-noise  time  delay  measurement.  Also  the  ellipses 
tend  to  align  with  the  tangent  to  the  measurement  lines 
defined  in  Section  IVrB  which,  in  this  case  of  time  delay 
measurements,  are  hyperbolas. 

2 .   Maneuvering  Target 

The  next  paragraphs  discuss  results  obtained  when  a 

selected   target  maneuver   is  simulated  ---  a  total  turn  of 

180  degrees  at  9  degrees  per  minute  with  a  constant  speed  of 
5.0  m/s. 

a.   One  Buoy  in  Action  at  a  Time  • 

Figure  ^1  shows  the  tracking  obtained  by  pro- 
cessing freguency  and  bearing  measurements  provided  by  a 
single  buoy  located  close  to  the  center  of  curvature  of  the 
target  path.  The  filter  takes  some  time  to  react  to  the 
target  maneuver  and,  when  it   does,   an   over-correction   is 
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applied.  The  result  is  that  the  estimated  path  is  very  dif- 
ferent from  the  true  path  of  the  target.  About  the  same 
result  is  obtained  when  the  buov  is  placed  on  the  other  side 
of  the  trackf  i.e.^  a  small  delay  to  react  to  the  maneuver 
followed  by  an  over-correction. 

In  Figure  ^2  the  freguency  measurements  were 
eliminated  and  only  bearing  indications  were  processed  to 
generate  an  estimated  track  very  close  to  the  one  of  Figure 
^1.  The  alignment  of  the  ellipses  with  the  measurement 
lines  is  again  very  clear  and  may  suggest  where  two  buoys 
should  be  placed  in  order  to  improve  the  tracking. 

b.   Two  Buoys  in  Action  at  a  Time 

Figure  ^3  shows  what  can  be  obtained  by  process- 
ing the  measurements  provided  by  two  buoys.  The  range  of  the 
buoys  was  adjusted  so  that  both  are  in  contact  with  the  tar- 
get during  all  the  path  shown.  Both  are  Difar  buoys  and  pro- 
vide bearing  and  frequency  measurements  with  a  period  of  100 
seconds^  which  are  processed  simultaneously  by  the  filter. 

In  Figure  ^^  the  measurements  were  processed 
separately  and  a  small  improvement  in  the  tracking  was 
obtained  at  the  end  of  the  pathr  although  during  the 
maneuver  the  simultaneous  processing  worked  better.  It  was 
observed  in  other  simulations  that  the  degree  of  improvement 
obtained  by  partitioning  the  bearing  and  the  time  delay 
measurements  is  dependent  on  the  maneuver  of  the  target   and 


la? 


B   100   100   S.0  d 
F   100   100   0*04h 


3.ee 


-4.ee     -3,38     -2.75     -3.13     -l.Sd 

(kn) 


-e.88     -C.eS     e,38        t.60 


.Eignre__4] jL_Irue.  .anri.. estimated _tra.j.eator.ies..  far  a. 

maneuvering  target  -  Simultaneous 
frequency  and  bearing  measurements. 


148 


3,0e 


B   100   100   5,0d 


a. 38 

1.75 

1.13 
(km) 

e.se 


-e.i3  « 


-0.75 


-1.38  _ 


-2.00 


<     t     t     1     t     I L- 

'4.00   -3.38  -2.75  -2.13   -1.50   -0.88   -0.25  0.38 

(km) 


1*00 


..Figure  42  ;  True  and  _  estimated .  traj.ec,toriea  .  f.or  .a 

maneuvering  target  -  Bearing  measurements 
only. 


149 


E.e0 

1,33 

0«75 

0.13 
(km) 

-0.50 
-1.13 

-1.75 

-2.3S 
-3.00 


B   100   100   5.0  d 
F   100   100   0.0^h 


i5«=i= 


w-^-&- "  ^^^..  ,'^^ 


1 


1 


1 


1 


i 


1 


-4.50   -3.88  -3.35  -3.63  -3.00   -1.38   -0.75  -0.13  0.50 

(km) 


F  i  gur e  43 


True  and  estimated  traj e^ctories  f or  a 
maneuvering  target  -  Simultaneous  pro 
cessing  of  frequency  and  bearing 
measurements . 


150 


the  relative  position  of  the  buoys.  The  worst  situation  is 
when  the  target  maneuvers  toward  the  buoys.  This  happens 
when  the  partial  estimates  are  worse  than  the  predicted 
values  and  occur  sometimes  during  fast  maneuvers^  or  slow 
maneuvers  with  long  intervals  between  consecutive  measure- 
ments • 

In  Figure  ^5  the  low-noise  time  delay  measure- 
ment was  added  and*  instead  of  reducing  the  errors*  the 
tracking  became  worse.  The  direction  of  the  ellipses  is 
clearly  along  the  tangent  to  the  hyperbolas  of  position* 
indicating  that  the  filter  is  uncertain  about  the  range  of 
the  target  to  the  buoys  during  all  the  maneuver.  This 
characteristic  of  the  measurement  lines  in  this  case*  and 
the  fact  that  the  target  is  modelled  by  the  filter  as  fol- 
lowing basically  a  straight  path*  explain  the  form  of  the 
estimated  track.  As  soon  as  the  target  reassumes  a  constant 
path  the  filter  starts  to  recover. 

In  Figures  46  and  47  the  iterative  techniques 
described  in  Section  IV*B*2*b  were  applied  to  the  frequency 
and  time  delay  measurements.  In  Figure  46  only  one  itera- 
tion was  executed;  in  Figure  47  the  iterated  formulas  were 
applied  three  times.  The  improvement  is  clearly  seen  by  com- 
paring with  Figure  45. 

As  with  the  partitioning  of  the  measurements*  it 
was  observed  from  other  similar  simulations  that  the  appli- 
cation of  the  iterative  techniques  does  not   always   provide 
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better  tracking.  It  depends  on  the  type  and  direction  of  the 
maneuvers*  the  position  of  the  buoys  and»  consequent  1 y »  on 
the  position  of  the  measurement  lines. 

In  the  case  of  the  non-maneuvering  target  the 
use  of  two  buoys  instead  of  only  one  was  almost  always  suf- 
ficient to  provide  good  tracking.  If  for  the  maneuvering 
target  one  uses  three  buoys  judiciously  positioned*  instead 
of  only  one  or  two*  one  can  also  obtain  good  tracking  as 
clearly  shown  in  Figure  48. 


C.   SUMMARY 

Below  are  listed  the  most  important  facts  observed  from 
the  simulations: 

---  The  tendency  of  the  error  ellipses  to  align  with  the 
measurement  lines  was  observed  in  Section  IV*8.  The  simula- 
tions show  that  the  sample  covariance  ellipses  also  follow 
this  t  rend. 

---  The  tendency  of  the  error  ellipses  to  align  with  the 
measurement  lines  can  be  of  great  help  in  the  practical 
solution  of  filtering  situations.  If#  for  example*  the  error 
ellipse  is  very  thin*  i.e.*  have  a  very  high  ratio  a/b*  then 
the  best  measurement  to  process  is  the  one  whose  measurement 
line  is  approximately  perpendicular  to  the  major  axis  of  the 
ellipse.  In  the  tracking  oroblem  this  can  be  obtained  by  a 
bearing  measurement  from  a  buoy  placed  along  the  perpendicu- 
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lar  to  the  major  axis  of  the  ellipse*  or  by  a  time  delay 
measurement  from  two  buoys  placed  along  the  direction  of  the 
major  axis  of  the  ellipse, 

---  The  frequency  measurement  have  a  very  small  effect 
on  the  tracking  accuracy. 

---  The  use  of  frequency  difference  or  frequency  ratio 
measurements  improve  computing  efficiency  (by  eliminating 
the  need  for  having  the  rest  frequency  as  a  state  in  the 
filter  model)  without  noticeable  effect  on  accuracy. 

---  The  partitioning  of  the  measurements  improves  com- 
puting efficiency*  may  provide  better  tracking  accuracy  and 
allows  the  processing  of  measurements  as  they  occur*  and 
thus  is  of  great  practical  importance. 

---  For  low  noise  measurements  the  iterative  techniques* 
although  requiring  an  increase  in  computing*  are  capable  of 
providing  better  tracking  for  maneuvering  targets, 

---  Tracking  with  only  one  buoy  is  acceptable  only  for 
non-maneuvering  targets, 

---  The  relative  position  of  the  buoys  is  one  of  the 
most  important  factors  which  influences  the  quality  of  the 
t  rack  1 ng. 
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VI.   CONCLUSIONS 


A.   SUMMARY  OF  RESULTS 

The  optimal  estimation  of  characteristics  and/or  parame- 
ters of  a  certain  class  of  nonlinear^  dynamic^  stochastic 
systems  has  been  studied  in  a  probabilistic  environment 
using  Bayes  formulation  concepts.  Approximate  solutions  and 
filtering  algorithms  were  generated  andf  among  them^  the 
known  Extended  Kalman  Filter  equations  and  higher  order 
filtering  equations  have  been  obtained  through  this  method. 

The  problem  of  tracking  submarine  targets  using  special 
passive  sonobuoys  was  modelled  and  with  this  model  extensive 
simulations  were  executed  to  allow  the  study  of  the  problem 
i  n  detai 1 • 

Most  of  the  results  indicate  that  the  frequency  measure- 
ments have  minimal  effect  on  the  filtering  process.  The 
small  contribution  to  range  information  that  thev  do  pro- 
vider when  associated  with  bearing  measurement s^  can  nor- 
mally be  obtained  otherwise  by  judicious  placement  of  subse- 
quent buoys. 

The  utilization  of  other  types  of  measurements  such  as 
the  frequency  difference*  the  frequency  ratio  and  the  time 
delay  measurement s*  proves  to  be  a  great  help  in  improving 
computing   efficiency  by  eliminating  the  rest  frequency  as  a 
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necessary  state#  and  thus  reducing  the  dimensionality  by 
one.  Alsof  especially  with  the  time  delay  measurement*  a 
great  improvement  in  tracking  accuracy  is  possible. 

The  concept  of  partitioning  the  measurements  and  pro- 
cessing them  separately*  even  if  they  occur  simultaneously* 
is  shown  to  bring  advantages  in  computing  efficiency  and 
also*  for  nonlinear  measurements*  in  tracking  accuracy. 
This  concept  is  of  great  practical  importance*  for*  in 
situations  where  the  measurements  are  sporadically  and  non- 
periodically  received  it  is  most  important  to  be  able  to 
process  the  measurements  as  they  naturally  occur. 

The  graphical  interpretation  of  the  action  of  Kalman 
filters*  developed  in  this  work*  provides  insight  into  the 
importance  of  each  variable  of  the  problem  in  the  filtering 
process.  The  direction  and  magnitude  of  the  correction  which 
is  applied  to  the  predicted  values  to  generate  the  new  esti- 
mate values  can  now  be  anticipated  as  a  function  of  the 
measurement • 

Nonlinearity  errors  have  been  graphically  presented  and 
iterative  techniques*  including  the  known  Iterated  Extended 
Kalman  Filter  equations*  have  been  suggested  to  counteract 
their  disruptive  effect.  The  application  of  these  tech- 
niques to  the  tracking  problem  shows  that  improvement  in 
tracking  accuracy  is  possible. 

The  graphical  interpretation  also  indicates  the  very 
practical   conclusion   that  the  error  ellipses  tend  to  align 
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with  the  measurement  Hnes»  as  defined  in  Chapter  IV.  This 
provides  guidance  for  optimal  positioning  of  the  buoys  and 
the  types  of  measurements  to  process.  This  observation  is 
reinforced  by  the  simulation  results  of  Chapter  V, 


B,   SUGGESTIONS  FOR  FUTURE  RESEARCH 

The  concept  of  partitioning  the  measurements  allows  one 
to.  process  the  measurements  separately  and  opens  some  ques- 
tions for  future  research^  such  as»  in  which  order  should 
nonlinear  measurements  be  processed  to  obtain  maximum  effi- 
ci  ency? 

The  graphical  interpretation  of  the  filtering  process 
allows  anticipation  of  the  direction  and  magnitude  of  the 
corrections  which  are  then  applied  to  predictions  to  gen- 
erate new  estimates.  This  situation  suggests  future  research 
in  the  determination  of  the  optimum  characteristics  of  the 
measurement  functions  which  in  turn  can  determine  the 
ODtimum  positions  and  characteristics  of  sensors. 

The  expansion  of  the  model  to  include  accelerations 
should  be  considered  if  maneuvering  targets  have  to  be 
tracked  efficiently  and  the  increase  in  computing  power  can 
be  obtained.  Consideration  of  the  depth  of  the  target  and 
the  inclusion  of  uncertainty  about  the  oosition  of  the  buoys 
are  additional  ways  to  extend  this  study. 
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The  model  and  simulation  implemented  herein^  without 
benefit  of  classified  information^  is  very  flexible^  but 
many  idealizing  assumptions  have  been  made.  Ihus/  a  natural 
extension  of  this  work  is  to  apoly  the  algorithms  and  con- 
cepts developed  to  a  real  problem  using  actual  sensor  data. 
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APPENDIX  A  :  PROBABILITY  RELATIONS 


The  rel at  i  on 


X  =  a  +  G. V 


(A.l) 


is  given  where  x  and  v  are  random  vectors  of  dimensions  n 
and  Qf  respect  i  vel  y  r  a^  i  s  a  n-di  mens  i  onal  deterministic  vec- 
tor and  G  is  a  (nxq)  matrix  of  deterministic  coefficients. 


The  joint  density  p  (v^)  is  given  and  the  joint   density 


p  ( X.)  ^  s  wanted. 


Case  ( i )  -If   n=q   andG    existsrthe   solution   is 
gi  ven  by  1 1 31 r 


h 


V    =    G  Tx    -    a]     =    f  (x) 


(A. 2) 


and 


P    (x) 

X  — 


=      det 


fci]  I . 


p     (f(x)) 
V 


(A. 3) 


Case  (ii)  -  If   n  <  a   one  could  add  dummy  variables 


n+1  n+1         n+2  n+2  q  q 


and    create 


x'    =    a'    +    G' .V 
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where 


Ui 


X 

L  qj 


a'  = 


Nowf  i  f  G  •"■'•  exists* 


=  C^  [x'  -  a'l  =  fix*) 


(A. a) 


From  Case  (i)»  the  solution  is 


P  ,(x')  = 

X 


det 


J-f 
3x'- 


.  p  (f  (x')) 

V  —  — 


(A. 5) 


and 


p  (x)  = 

X  ~ 


p  ,(x').dx  ,-.dx  ,„,..dx 
X  —      n+1    n+2      q 


(A. 6) 


Case  ( i  i  i  )  -  If    n   >   q    the   following   development 
appl i  es. 

Let 


16a 


X       = 


a    = 


q+1 


q+1 


q+1 


G    = 


^qti.i W}'^ 

^,1  n,q 


^+1 


From    (A.l),    one    also    has    that 


x'    =    a'    +    G' .V 


and    then*     if    G*         exists* 


V    =    G  '"^    [  X  •    -    a  •  1     =    f  (  X  •  ) 


(A. 7) 


From  Case  (i)  one  then  has 


P  ,(X,/X.f..,,X   ) 

X  ■  1   2      q 


=   det 


8x'- 


.  p  (f(xM) 

V  —  — 


(A. 8) 


but 


Now  consider  the  variable  x   .  Given  (A, 7)* 

q-lfl 


P(X^-|X,X»...,X       )      =      D(x!VfVf...fV      ) 

q+1      12  q  q+1      12  q 


q+1  q+1         q+1  J.  1  q+1, 2    2  q+1  ,q    q 


=    a         +    g  .V 

q+1        ^+1     - 


or 
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^q+1     =    ^q+1      ^    Vl 


.f(x')     =     1 


q+1 


(A. 9) 


From    (A. 9)    comes 


q+1        12  q 


''\.i  -  Vi' 


and 


p(x     fX     /...fX  )     =       6 (x  -    1  ).p     (x     fX     ,...,x     ) 

1        2  q+1  q+1  q+1  x'      1       2  q 


Doing  the  same  thing  for  x  ^,     »,,,    x    #   one   finally 

q+2  n 


gets 


p    (x)    = 

X  — 


6(x.     -    1 . )     .       det 
1  1 


fo] 


.    p    (f(xM) 
V    —   — 


i=q+l 


(A. 10) 


For  the  special  case  of  q  =  1/  the  vectors  _x  '  ^   a^'   and 
g'  become  scalars  and  the  sbove  equations  lead  to 


f (x')  =  f (x  )  =  (x   -  a  )  /  g 

1       1     11. 


tl-  '  '  ^1      ' 


1   =  a   +  (x   -  a  )g  /g 

r     r      1     1   r   1 


and  Equation  (A, 10)  becomes 
n  , 


p  (x)  = 
X  — 


6(x.  -  1.  )  . 

X       1 


1/g 


•p((x   -a)/g) 

V    1      11 


1=2 


(A. 11) 


Another  similar  formula  can  be  obtained  for  this  soe- 
cial  q  =  1  situation.  The  matrix  G  in  (A.l)  is  now  a  n- 
dimensional  vector  and^  if   one   breaks   (A.l)   into   its   n 
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individual  scalar  equations  one  can  obtain 


p(x^lv)  =   6(x^  "  ^i  '  ^i^' 


Also    note    that 


p(x^; X     , v) 


p(x^  I  v) 


i    P    J 


and 


p(  X  .  »  X    .  I  V ) 

1     J 


p(x.|x,rV).p(x      Iv) 

1       J  3 


p(x     !v).p(x     Iv) 

1  J 


And    from    these    relations* 
n .  n 


p    ,  (xlv)    = 

x/v 


p(xlv)     = 


6  (  X       -    a       -    g    V  ) 
i  i  i 


i=l 


i=l 


and    f  i  nal 1 y # 


p    (x)    = 

X 


6(x.    -    a,    -    g.v).dv 
1  11 


(A. 12) 


i=l 
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APPENDIX  B  :  COMPONENTS  OF  COVARIANCE  MATRIX 


For  a  single  measurefnent  of  the  form  z  =  H.x  ♦   v#   the 
observation  matrix  becomes  the  vector 


H—  th    h„   •••   hi 
-  1    2         n 


where  n  is  the  dimension  of  the  system. 


The  Kalman  gain  is  given  by 


T      T     2 
G  =  P'H  [HP'H   +  r  J 


(B.l) 


where  P*  is  used  to  represent  the  prediction  error  covari- 
ance  matrix,  which  is  symmetric,  and  r  is  the  standard  devi- 
ation of  the  measurement  error. 

T 
Since  H  is  in  this  case   a   vector,   the   product   HP'H 

will  be  a  scalar. 


Let 


C  r  [HP'H^  +  r^f-^ 


(B.2) 


The  vector  P'H   can  be  shown  to  be 


P'H   = 


n 


,2,  p.'.  h. 

n 

I      P'.  h. 


Z   P'.  h. 


and  thus  the  constant  c  is  given  by 


(B.3) 
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n   n  2 

C  =  1  /  I  Z   Z  P'  h  h    +   r^l 
i=l  j=l  ij   i  j 


(8.4) 


From  (B.l)  one  now  sees  that  the  Kalman  gain  is  a  vec- 
tor whose  comoonents  are  the  components  of  P ' H^  multiplied 
by  the  constant  Cf  i.e.  the  components  of  the  gain  vector 
are 


n  n   n  o 

g    =  .Z-   p'  h.  /  I  .Z,  .Z,  P.'.h.h.  +  r^] 
g     3=1   mj  J      1=1  J=l  ij  1  J 


(B.5) 


The  estimation  error  covariance  matrix  is  computed  from 
the  equation 


P  =  II  -  G.H]  .P'  =  P'  ♦  AP 


(B.6) 


where 


aP*  =  -  GHP'  =  -  cP'H  HP' 


(B.7) 


Let • s  def  i  ne 


n 
a.  =  .Z,   P*.  .h. 
1   J=l    ij  3 


=  ^n^  '  ^iz'z   *   '• 


+  p'  h 
in  n 


(B.8) 


Then»  f  rom  (B. 3) ^ 


HP*  =  tP'H^J  ^  =  [a,   a-,   ...   a  1 

12         n 


(B.9) 


and  f  rom  (B. 4 ) , 


n 


c  =  t.E,  h.  a  ,  +  r^J 
x=l   X  1 


2.-1 


(8,10) 


From  (B.7)  the  increment  aP'  to  the   error   covariance 
mat  r i  X  is  then 
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A  P'=  -  c 


a  a 

1   2 


a  a 

1  2 


a.  a 

1  n 


a_a 
2  n 


a  a 

1  3 

a  a 

2  3 


a„a 

3  n 


a  a 

1  n 

a  a 

2  n 


"J 


(B.ll) 


The  individual  terms  of  the  estimation  error  covariance 
matrix  wi]1  be  given  by 


p    =0    =  p'   "  c.a  .a   =  D     " 
ij     ji      ij        i   j     ij 

n  n  o  ^ 

-a, a   /(lEp'hh    +   r^) 
i   j     i=lj=l  ij   i  j 


or 


13 


n  n 

k=lnSl 


ilnSl  ^^^Pij^' 


'^-;^' 


0.'.  r^ 


n   n 
k=l  ni=l  k  m  k 


+   r 


(B.12) 
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APPENDIX  C  :  CHARACTERISTICS  OF  MEASUREMENT  FUNCTIONS 


in  Figure  C-1  some  curves  of  the  form  f(x) 


constant 


are  shown^  where  f  ( x.)  is  assumed  given.  The  vector  2L  arid  a 
constant  value  d  are  also  given.  The  vector  x^  is  defined 
by  the  equation 


b  .   a  ^ 


(C-1) 


or 


x^  +  h.  (d/lhl^l 


(C-2) 


where  h^  is  the  gradient  of  f  ( £)  with  respect  to  x_   taken   at 
point  (x  tx    ),    and  h/!h|  is  a  unit  vector  along  h. 


_h,  =  h.  e^,  ♦  h  ^  e 


2-2 


4"-'^ 


Ihl  =  (h^  ♦  h2jl/2 


•^2  =4"-'' 


lo  determine  the  location  of  j<  numerical  data  is 
required.  However,  it  is  possible  to  determine  qualitatively 
whether  x.^  ends  on  the  curve  fix^)  +  d»  or  on  one  side  or 
the  other  of  this  curve.  This  is  done  in  the  following 
paragraphs  for  three  special  functions. 

Case  ( i  )  -  For  linear  f  ix^) ,  j("  always  ends  on  the  curve 
fix)    =flx^)  +  d,  as  shown  below 
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Figure  C  -  1  :  Problem  geometry 
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Let  fijt_)    =  px   +  qx 


I  hen  h   =  p 


»    h   =  q   and  jhp  =  p^  f  q^ 


From  Equation  (C-2)f 


X       '    }^     *    (pe   +  qe  ).d/(p^  +  q^) 


so 


x^   =    x^   t    pd/(p2    +    q2) 
X2   =    x^   +    qd/(p2    +    q2) 


NOM 


f(^^)    =    P^     +    qx^     =    px^     +    p2q/(p2     t    q2)    + 


♦    qx^     f    q^d/Cp^     +    q2  ) 


=    px^     *    ^'^l     *    cl    =    f(x_^)     +    d 


Case  ( >  i )  -  For  fix)    defined  as  belowf  x°   always   ends 
on  top  of  f  ( x^  )  f  df  as  in  the  Case  (i). 


fU)    =  [(x^  -  p)2  +  (x   -  q)2)l/2 


Now 


h,  =  (x,  -  p)/m 


h   =  (x^  -  q)/m 
2      2 


where 


m    = 


i^^l 


p)2  t  (x^ 


-  q)2jl/2=  fCx^) 


and 


•h'2 


h2  t  h? 


=  1 
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Hrom  (C-2), 


X.   =  X.   "•"  il'd  =  J<   ♦  l(x   -  P)e   +  ( x^  -  q)e_].cl/m 


so 


b  _ 
b  _ 


X   +  h.  d  =  X-   +  ( 


\ 


-  p)d/m 

-  q)d/m 


Now 


fCx^J 


[(x^  ♦  h  d  -  p)2  +  ()^   ♦  h  d  -  q)2]^2_ 


=  [(x^  -  p)2  f  2h  d(x^  -  p)  f  h^d^t 

+  (x^  -  q)2  +  2h  d(x^  -  q)  +  h^  d^]^^^ 
2  2    2  2 


Using  the  values  of  h   and  h  » 

f(x_^)  =  [(x^  -  D)2  +  (x^  "    qf-     +  d^  + 

+  2d[(x^  -  p)  /m  +  (x^  -  q)  /m]J^/2; 
1  2 

=  [f(x^)^  +  2d.f(x^)  +  d^]-^/^= 


=  [(f(x^)  +  d]?  J^^^  = 


=  fCx**)  +  d 


tase  ( i  ii )  -  For  f(x)  as  definea   below   and   shown   in 

b  a 

Figure   C-2#   x   always  ends  before  the  curve  f(x  )  +  d/  for 


(d|       <     TT  . 


f(x_)    =    arctan[(x       -    q)/(x       -    p)] 


17a 


(x)  =  f  (x^)  +  d 


Figure  C  -  2  :  Geometry  for  Case  (iii) 
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The  case  j  d  I  >  tt/2  is  not  of  interest  since  the  lines  h 
and  f  I  x_J  =  i  (.j^  )     +  d  are  diverging.  For  jdj  <  tt  /2r 


-  (x^  -  q)/n 


(x^  -  p)/n 
1 


where 


n  =  (x^  -  p)2  +  (x^  -  q)' 


and 


1^1 


2  _ 


h2  ♦  h2 


=  1/n 


h rom  Equation  (C-2)f 


x^    =    x^    f    d.  (-(x^     -    q)e       +    (x^ 

-  -  2^-1  1 


-    p)e    J 

-2 


SOf 


11     "         2    "    ^ 


Now 


f(x^J    =    arctanKx^    -    q)/(x^    -    p)]     = 
-     -  2  1 


arctan     [(x^    -    q)     +     ( x^    -    p)d]/ 
2  1 


l(x^     -    D)     -    (x^     -    q)dl 
1  2 


or 


f(x°J    =    arctanKs    +    d)/(l    -    sd)] 


where 


s  =  Ix   -  q)/(x^  -  p) 


=  tan  [f (x^)] 
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SOf 


tan  [f (x  ))  =  (s  +  d)/(l  -  sd)  = 


(tanlf(x^)J  +  dJ  /  tl  -  d.tan  (f (x^)J 1 


=  tantfCx  )  +  arctan(d)) 


OP 


f(x^J  =  f(x^)  +  arctan(d) 


Since  |d)  <  i\  /2,    then  arctan(d)  <  d   and  tt       will   not 


reach  the  line  f(x  )  f  d. 
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APPtNUlX  D  :  ANALYTICAL  EVALUATION  OF  KALMAN  FILTERS 


After  the  design  of  a  filter  an  analysis  phase  is 
necessary  to  verify  its  behavior  in  various  representative 
si  tuat  i  ons • 

A  "filtering  situation"  is  considered  here  to  be  com- 
pletely specified  by: 

---  the  true  parameters  and  initial  conditions  of  the 
system  which  generate  a  unique  track  x(0)r  x(l)r  ..,/  x(t  ), 

---  the  true  parameters  of  the  sensors  and  their  meas- 
urement schedules, 

---  the  parameters  and  initial  conditions  assumed  by 
the  filter  for  the  system  and  for  the  sensors. 

Ihe  approach  normally  used  to  determine  filter  behavior 
is  to  simulate  the  desired  situation  and  execute  hundreds  or 
thousands  of  Monte  Carlo  runs  and  compute  sample  statistics. 
The  most  useful  results  of  this  process  are: 


m 


a^  IkJ   -   the  sample  mean  of  the  estimated  state  vector 
at  time  t,   obtained  from  m  Monte  Carlo  runs. 


^°lk)  =  _a™(k)   -  j<(k) 


the   sample   mean   of   the 


estimated  error  vector  at  time  t,  . 

k 


m 


b  Ikj   -   the  sample  second  central  moment  of  the  state 


m 


(or  errorj  estimates  about  a  (k)^  at  t,  . 

...      —  k 
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Ihe  objective  of  this  section  is  to  study  the  possibil- 
ity of  obtaining  values  of  £(  k ) ,  £(k)  and  b^(k)/  which  are 
approximately  the  values  of  ^  (.k)  ,  £™(k)  and  b"^  ( k )  when  m  is 
very  larger  without  the  use  of  the  time  consuming  Monte 
Carlo  operat  i  ons  . 

1  ---  Theoretical  Solution 


Unce  a   "filtering   situation^r   as   stated   abover  is 

defined/   the   operation   of  a  Kalman  Filter  is  described  by 

the  recursive  equations  below  where  only  the   dependence  of 
each  term  on  the  previous  estimates  is  stressed. 


i(kf  Ilk)  =  f (x(k)) 


(D.l) 


H(ktllk)  =  ())(i(k)).P(k).i|)Mi(k))  + 


t  g(i(k)).Q(k).g^{i(k)) 


(0.2) 


G(x(k+1 |k),P(k+l ik))  =  P(k+1 Ik) .H^ (x(k+l ik) ). 


lH(i(kf  1  |k)).P(k+l  |k).H'^(i(k+l  Ik))  +  R(k)]  ^ 


(D.3) 


x(k+l)  =  x(k+l|k)  ♦  G(x(k  +  1 lk),P(k+l Ik)).  (h(x(k  +  l))  + 


+  v(k  +  l)  -  h(x(kf 1  Ik))) 


(D.a) 
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P(ktn  =  (I  -  G(x(k+1 !k),P(k+l Ik)) .H(x(k+1 |k) )J .P(k+1 ;k) 


(0.5) 


In  this  case  these  equations  describe  a  deterministic 
process  were  it  not  for  the  random  measurement  noise  v(k+l)» 
the  only  random  input  to  the  "filter  dynamics"  since  a 
unique  track  is  considered.  The  initial  condition  of  the 
filter*  x^lO)  and  P(0)  are  also  determi  ni  st  i  cal  1  y  given. 

tquations  (0.^)  and  (0.5)  can  be  rewritten  as 
x^(k  +  l)  =  (J)  (x^(k  +  l  |k),P(k  +  l  Ik))  +  'I'  (x_(k  +  l  |k),P(kf  1  Ik))  .v(ktl) 

P(k  +  1)  =   f  (x^(k  +  l|k),P(k  +  l  Ik)) 

where  <j)  r  "F  and   r  are  generally  nonlinear  matrix   functions 
of  the  variables  indicated. 

Considering  (0.1)  and  (0.2)  one  can  then  write 


x.(kfl)  =  f_*  (£(!<), P(k))  +  g^*(^(k),P(k)).v_(kfl)    (0.6) 


P(k+1)  =  h*(x(k),P(k)) 


(0.7) 


JL  JL  JL 

'where  L  '  9.  '  i!  ^^^  also  generally  nonlinear  matrix  func- 
tions. Prom  (0.6)  and  (0.7)  one  can  clearly  see  now  thatr 
i  f  v^  i  s  a  discrete  white  Gaussian  noise  process*  the  joint 
process  <_x/P>  is  Markov, 

If  one  considers  a  new  vector  y^  formed  with  the  n  ele- 
ments  of  X  and  the  n,(ntl)/2  distinct  elements  of  the  (nxn) 
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symmetric  matrix  P,  y  is  an  tn.  (n  +  3) /21 -di mensi onal  process. 
Equations  (D.6)  and  (D.7)  can  then  be  combined  into 


y(ktl)  =  A(y(k))  f  B ( y ( k ) ) . v ( k+ 1 ) 


(D.8) 


where  A  is  a  vector  function  and  B  a  R   to  R™   matrix   func- 
tion/ m  =  n  ,  (n  +  3) /2, 

fhe  complete  behavior  of  the  filter  can  be  described  by 
the  probability  law  of  ^t    which  is  obtained  from 


p(^(k+l))  = 


p(^(ktl)  ;_y(k)).p(y(k)).d_y(k)      (D.9) 


and  p  (y^tkt  1 )  j_y  (  k  )  )  is  obtained  from  p(v(k  +  l))  and  Equation 
(D.ti)  in  the  manner  shown  in  Appendix  A,  The  initial  value 
^(OJ  is  deterministic  since  it  contains  the  initial  condi- 
tion of  the  filter  which  is  given, 

i.    ---  The  Linear  Case 

Kor  the  special  case  of  linear  dynamics  and  observa- 
tionSf  the  solution  can  be  obtained  in  a  simpler  way.  In 
this  case  the  functions  (j)  '  Q  and  H  ^r^  not  functions  of  the 
state  estimates  and  so  neither  is  G, 

Equations  (0,1)  -  (D,5)  can  be  grouped  now  into 

£(k  +  l)  =  <j)  (x^(k  +  l|k),P(k  +  l,k))  +  ^  (P(k  +  1  |k)).v^(k  +  l) 

PCk+l)  =  r  (P(k+1 !k)) 

and  finally/  in  a  recursive  way/ 
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x(k+l)  =  f *(i(k),P(k))  +  g*(P(k)) .v(k) 


(D.IO) 


P(k+1)  =  h  (P(k)) 


(O.ll) 


From  (D.ll)  it  is  clear  that  the  estimation  error 
covariance  matrix  is  now  a  deterministic  process  and  can  be 
precomputed*  together  with  the  gains. 

Equation  (D.IO)  then  becomes 


i(k  +  l)  =  t*(i(k))  t  a*(k  +  U.^(k+n 


From  (D.l)  -  (0,5)  the  values  of  f_      and  g*  can  be  found 
to  be»  for  this  linear  caser 

i(k  +  l)  =  II  ♦  G(k  +  n.H(k+l)J  ,<|)  (k).^(k)  + 

f  G(k  +  l).H(k+l).x-(k+l)  f  G(k  +  l).v(k+l) 


or 


x(k+l)  =  S(ktl).x(k)  +  F(k+1)  ♦  G(k+l).v(k+l)    (D.12) 


where 


S(k)  =  [1  +  G(k).H(k)]  .(})(k) 


F(k)  =  G(k).H(k).x(k) 


I  he   sought   after   behavior   of   the   filter   can    be 
described  by  the  moments  below/  which  agree  with  1121 
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a(k  +  l)  =  Eli(ktl)]  =  S(k  +  l)  .E  [x(k)J  t  F(k+1)  f 


+  G(k+l).E[v(k+l)] 


(D.13) 


e(k-H)  =  a(ktl)  -  x(k  +  l)  =  S  (  k  + 1 )  .£  li  (k)  1  ♦ 


+  II  -  G(k+l).M(k+l))  .x(k+l)  + 


t  G(k+l).Elv(k+l)]  = 


=  S(k+l).e(k)  f  D(k+l).x(k+l)  -  S(ktl).x(k)  ♦ 


tG(k+l).E[v(k+l)] 


(D.ia) 


b(k+l)  =  S(k+l).b(k).S^(k+l)  + 

+  G(ktl),Etv(ktU.v^(k  +  l)]GT(k+l)    (D.15) 

S    ---  General  Case  ( 

For  the  general  nonlinear  problem  one  returns  to  Equa- 
tion (D,9j  which  cannot  normally  be  solved  in  closed  form. 

Numerical  techniques  and  approximations  would  be  neces- 
sary that  may  use  more  computing  power  ana  present  worse 
results  than  the  Monte*Carlo  process  that  we  are  trying  to 
avoi  d. 

One  way  to  simplify  the  problem  is  to  assume  that   y(k) 

has   an  approximate  Gaussian  di st r i but i on#  even  with  all  the 

non I i near i t i es  shown.  Linearizing  Equation  (D,8)  around   the 
mean  value  of  y(k)  would  give 
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^(k  +  l)  =  A(y(k))  ♦  ^  A(jy(k)).  [y(k)  -  y(k))   + 

dX  —        — 


+  B(y(k)).v(k+1) 


(D.16) 


Taking  ^(0)  as  the  determi ni st i ca  11  y  known  initial 
value  y^iO)f  and  applying  the  expectation  and  variance  opera- 
tors to  Equation  (0.16)^  one  gets  the  recursive  equations 
for  the  moments  of  y(k): 


y(k  +  l)     =    A(y(k))     t    B ( y ( k ) ) .E  (v ( k+ 1 ) 1 


(D.17) 


Varly(k  +  l)J     =     I  ^T"  A  (y  (  k  )  ) )  .  Var  ly  (  k  ) )  .  [  -i-     Myi\c))f     ♦ 

3X       -  -  " 


3^ 


+    B(y(k))  .Var  [v(k  +  l)J  .BMy(k)) 


(D.18) 


Ihe  moments  we  are  interested  in,  a^r  e_  and  'b  are 
directly  obtained  from  subvectors  and  submatrices  of  the 
above  moments. 

The  primary  difficulty  with  this  approachr  howeverr  is 
in  obtaining  the  functions  A  and  B,  This  can  be  seen  for  the 
very  simple  case  belowr  as  an    example. 

Suppose  a  scalar  system  with  a  single  observation  is 
described  by  the  equations 

x(k  +  n  =  x^(k)  +  w(k) 

z(k)  =  x(k)  ♦  v(k) 
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tquations  (D.l)  -  (D.5)  give 


x(k  +  l  Ik)  = 


i^Ck) 


(|)(k)  =  2x(k) 


p(ktllk)  = 


Hi  ^(k).p(k)  f  q^(k) 


U(ktl)  = 


Hx^  (k).p(k)  +  q^(k) 
ai^(k).p(k)  +  q^(k)  +  r^(k) 


and/  after  the  appropriate  steps* 


x(k+l) 


p(k+l) 


i^(k)  + 


ax^(k).p(k)  +  q^(k) 
ax^(k).p(k)  +  q^(k)  +  r^(k) 


.  [x(k)  -  x(k)3 


iai^(k).p(k)  +  q^(k))  .r^(k) 
ai^(k).p(k)  +  q^(k)  +  r^  Ck) 


ai2(k).p(k)  ♦  q^(k) 
qx^(k).p(k)  f    q^(k)  +  r^(k) 


.v(k+l) 


The  mean  value  given  by  Equation  (D.17)  can  be  obtained 
in  a  simpler  way/  however.  It  is  the  result  obtained  by 
running  the  filter  in  the  simulated  environment  but  making 
the  measurement  noise  assume  its  mean  value  Etv(k)]/  nor- 
mally zero. 

If  one  uses  more  terms  in  the  expansion  of  Equation 
(D.8)/  better  results  will  *  be  obtained  at  the  expense  of 
increased  computing  effort  and  time. 
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Fop  this  general  case  it  is  suggested  that  these  equa- 
tions be  applied  to  a  simple  scalar  or  two-dimensional  non- 
linear problem  and  that  the  computing  power  required  to  find 
the  sought  after  moments  be  compared  with  that  required  to 
run  a  Monte  Carlo  process  yielding  the  same  level  of  accu- 
racy in  the  results. 
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APPENDIX  E  :  COMPUTER  PROGRAM 


Ihe  basic  structure  of  the  computer  program  written  for 
the  simulation  of  the  tracking  problem^  and  the  evaluation 
of  models  and  filtering  algorithms/  is  shown  in  Figure  E-1, 
The  solid  lines  represent  the  normal  flow  of  control  within 
the  program;  the  dashed  lines  show  the  extra  pathes  that 
become  available  whenever  an  interruption  is  requested  by 
the  user. 

After  a  brief  introduction  to  the  program  a  table 
called  MtNU  is  presented  to  the  user,  as  shown  in  Figure  E- 
2.  If  the  user  chooses  actions  number  1  or  5  the  result  is 
c 1  ear • 

Choice  number  3  provides  access  to  all  the  problem 
variables.  Since  the  number  of  variables  that  characterize 
each  simulation  is  very  larger  a  simple  way  to  deal  with 
these  variables  had  to  be  devised.  As  it  is  imolemented/ 
the  program  always  "remember"  the  values  given  to  the  vari- 
ables at  the  last  time  the  program  was  used»  so  that  only 
the  variables  to  have  their  values  modified  have  to  be 
addressed.  The  change  of  value  of  a  variable  is  simply  made 
by  selecting  the  page  where  it  appearsr  typing  the  letter 
associated  with  thevariable  and  the  new  desired  value.  The 
program  immediately  responds  by  presenting  back  the  entered 
value.    A   more   detailed   diagram  of  the  actions  resulting 
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I 


Tut  or  i  a  1 

In  format  i  on 


M  p  N  II 


Chanae 
P  rogram 

Cha  rac  t 
CPC 


r  hanae  Va \ ue 
o  f  Prob 1  em 
Variables 


STOP 


Simulation 


Printing  and 
^1 ot  t  i  nq 
Pout  i  nes 


Output 

f^  i  1  es 


F  i  au  re  E -  1 


Basic  s^^uct•u^e  of  oroara-T' 
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MFMU 


Press  the  nufrber  corresoonding  to  the  action 
des  i  red  and  c/r . 


I   -   PRFSENT  TUTORIAL  IMFORMATION 


2   -   ^^OniFY  PROGRAM  FLAGS 


3   -   FORMULATF  OR  MODTFY  THE  PRQBLtW 


a       -   START  THF  SOIUTIOM  0^  THF  PROBLE.M 


5   -   FNH  THE  PROGRAM  AMD  «^XTT 


Fioure    E-<£     ;     The    MP(>iU    tabi 


189 


from  choice  number  3  is  oresented  in  Figure  E-3.   A   typical 
page  of  variables  is  shown  in  Figure  E-^. 

Choice  number  ^  simulates  the  tracking  problem  accord- 
ing to  the  values  defined  in  the  previous  step,  A  block 
diagram  of  the  actions  involved  in  the  simulation  is 
presented  in  Figures  E-S*  E-6  and  E-7, 

Initially  all  the  important  problem  variables  are 
printed  with  their  current  values  in  a  form  as  shown  in  Fig- 
ures E-tt/  E-9  and  E-10,  During  the  first  run  the  scheduling 
of  all  the  events  involved  in  the  simulation  is  also  printed 
for  future  reference^  as  shown  in  Figure  E-11. 

With  choice  number  2  from   the   MENU   a  new   table   is 

presented   as   shown   in   Figure   E-12.   This  table  in  also 

presented  to  the  user  whenever  he  requests  an  i nterrupt i on^ 
at  any  time  or  point  within  the  program. 

The  extra  choices  that  now  become  available  are  almost 
self-explanatory  but  it  should  be  added  that  with  choice 
number  7  the  simulation  is  ended  and  the  partial  statistics 
are  computed  and  written  in  the  appropriate  output  files? 
choice  number  6  allows  the  modification  of  any  variable  in 
the  middle  of  the  simulation^  in  the  same  way  as  explained 
before?  choice  number  0  transfer  control  back  to  MENU  or  to 
the  point  of  interruption. 

A  listing  of  the  program^  in  C  language/  is  available 
at  the  Electrical  Engineering  Department/  Naval  Postgraduate 
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Sel ect 
Variable 
To  Be 
f^od  i  ^  i  ed 


Chance 
Value  of 
Se 1 ec  t  ed 
Variable 


*  ronn 
CPC 


from 


P  resent 

P"  i  rsf  Paqe 

of  Variables 


Choose   Action 


Print  Pane 
If  It  Was 
Modified 


to 

CPC 


Print  P  aqe 
If  It  /las 
^od  i  f  i  ed 


Present 
^fext  Paqe 
of  Variables 


to 

MEMU 


Fioure  E-i  :  Chaoqino  ^he  value  of  problem  variables 
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To  change  the  value  of  any  variable^  type  the  indicated 
letter^  at  least  c^e  soacer  the  desired  value  and    c/r. 
To  see  the  other  variables  of  the  nroblenn,  tvpe  1  c/r. 
To  go  back  to  menur  tyoe  0  c/r, 

(  1  meter  =  1.09^  yards;   1  rreter/sec  =  K^^q  knots  ) 


TARGFT  INITIAL  PARAMFTFRS 


initial  x    -  nos  i  t  i  on 
y  -  Dosition 
speed 
headi  no 
f  reauenc  v 


0.0  km 

0.0  krr- 

7,0  m/sec 

30S.0  deq 

SOO.ft  hert  z 


a 
b 
c 
d 
e 


standard  deviation  of  ^orcino  functions 


speed/sec 
headi  no/sec 
f  reauenc  v  /  «;ec 


0,0  m/sec/s«c  f 
O.n  deg/sec  q 
<}  ,0       hert  z/sec     h 


Fiaure  E-^  :  Typical  oaae  of  variables 
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from  MFNU 


Print  Value  of 
All  Variables 


T 


run  =  t 


1 


Initialize  Variables 

and  Random  Number  Gen 
for  a  Mew  Pun 


t  i  rre    -     I 


Is  taraet  to  be  uodated 
every  second? 


Are  There  Anv 
Other  Fvent  s 
Occu  r  r  i  no  a  t 
This  Time?  ' 


Print  Table 

of  Events 


I UoHate 
^  Target 

Parame t  e  r s 


Figure  E-b 


Simulation  T 
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t  i  Te  =>  1 


Is  time  >  duration 
of  a  run? 


Is     run    >    nurrber 
of     runs? 


run  =+  1 


CoTDu  t  e 
Statistics 


Senci  Results! 
To  0 u  t  D u  t 
^^  i  1  e  s 


to 
1.  MENU 


Fioure  E-b  :  Simulation  TI 
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Chanae    Taraet 
Pa  rannet  ers    to 
S  i  Tiu  1  at  e 
Maneu ve  r 


Is    t  aroet     to    be 
upfiaterf    before 
•measurement  s?  — ' 


/' 


Oe t  e  rm  i  ne 
Pred  i  c  t  ed 
State 


Se 1 er  t  Ac  t  i  on 
FroT  Table 
Of  Events 


Is  the  reaues  t  ed 
number  of  iterations 
satisfied?    • 


Simulate 
Sel er ted 
Measurement 


Process 
^easu  remen  t 
To  Oh t  ai  n 
fs*"  i  "^at  e 


Store  Values 
of  Paramet  ers 
and  Fs t  i  mat  es 
for  Statistical 
Purposes 


Target 

Pa  rame t  ers 


i  Linearize 
.^.  About 

Predi  c  t ed 
State 


1 nea  r 1 ze 
'^bout  Part  i  a  1 


I  Fs  ♦•  i  ma  t  e 


F  i  au re  E - 7 


Simulation  Til 
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TARGFT     IMITIflL    PAhfAMFTFRS 


initial 


X  -  DOS  1 t ^ on 
V  -  oosition 
speeH 
headi  nq 
f  requenc  v 


0,0  Km 

0.0  km 

7,0  m/sec 

"^OS.O  dpq 

'^OO,^  hertz 


a 
b 

c 
d 

e 


standard  deviation  of  forcina  functions 


soeed/sec 

headino/ser 

fr«?quencv/sec 


0,0  m/sec /sec  f 
0.0  deq/sec  c 
0.0   hertz/sec     h 


GENERAL  niRtCTIONS 


nuT^ber  of  HIFAR  buoys, 
number  of  LOFAR  buoys 

number  of  maneuvers 
number  of  runs 
durat  i  on  of  runs 
number  of  key  ooints 

seed 

averaqe  sound  velocity 

range  of  the  buoys 


^ 

buoys 

a 

0 

buoys 

b 

? 

man 

c 

?0 

runs 

d 

aooo 

sec 

e 

^0 

po  i  n t  s 

f 

1^3^^ 

a 

1500,0 

m/sec 

h 

8.0 

km 

i 

FILTER  PIITIAL  PARA^lFTFRS 


initial  x  -  oositinn 
y  -  nosition 
^oee'l 
^  e  =»  a  i  n  n 
f  r  e  n  1 1  e  n  c  V 


-0,S  km 

n  .0  km 

6.0  T  /s'='c 

^lO.n  n-;. 

S  0  0 , S  hertz 


a 
b 

c 
d 
e 


standard  deviafion  of  forcina  functions 

spp»ed/sec  :  0,00S  m/sec/sec  f 

heaainn/ser  t  0.050  aeq/sec  n 

frequency/sec  '     0,0  10  hertz/sec  h 


Fiaure  fc-b  :  List'ina  of  oroble"'  variables  T 
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TNTTTAl.     COVAPIfli\CE    MATRIX 


X  DOS 


V  P05 


SDPer! 


heaoo 


f  reo 


X  pos 

l.O 

a 

0.3  b 

0.0 

c 

0.0 

(1 

.  0.0 

e 

y  POS 

* 

I  .U  q 

0.0 

h 

o.u 

i 

0.0 

J 

speed 

* 

* 

2.0 

m 

0.0 

n 

0.0 

o 

heada 

* 

A 

* 

100.0 

s 

0,0 

t 

f  req 

* 

* 

* 

A 

1.0 

y 

(km,  m/sec»  Hearees  and  hertz  crossmulMolied) 


MtASURFiMFUTS  SCHFDMLF   I 


huoy 

t  yoe 

start 

period 

mean 

1 

B 

100 

1  un 

0.0 

1 

F 

mo 

100 

0.0 

? 

8 

500 

ion 

0.0 

? 

F 

500 

1  0  0 

0.0 

s  t  dev 

5.0  abcaef 

0  .Oa  Qhi  j  k 1 

5.0  mnoDor 

O.Oa  stuvwx 


(start  and  period  in  secondsl 

(mean  and  std  dev  in  dpqrees*  hert?  or  seconds) 

(period  rrust  Oe  ?ePO  i^  Teasurement  is  not  used) 


MEASURFMFfJTS  SCHEDULF   II 


buoy 

t  vpe 

start 

pe  r  i  od 

mean 

1? 

T 

500 

1  00 

0.0 

13 

T 

500 

1  00 

0.0 

23 

T 

50  0 

1  00 

0.0 

3 

B 

SOO 

1  on 

o.n 

st  dev 

0.01  abcdef 

0.01  nhi  j  H 

0.01  nnopor 

5.0  stuvwx 


(start  and  period  in  seconds) 

(mean  and  std  dev  in  deqreeSf  hertz  or  seconds) 

(period  must  oe  zero  i^  measurement  is  not  used) 


F  i  oure  E -9 


Listina  o^  problem  variables  II 
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1st 


TARGET 

yiAfNjFuvFP 

5 

maneuver     : 

time     - 

60  0 

sec 

a 

vdot     = 

0.0 

n>/sec/mi  n 

b 

hdot     = 

-9.0 

d  e  q  /  m  i  n 

c 

frJot     = 

0.0 

hert  z /m  i  n 

d 

?nd 


TARGFT 

MAfsiFUVERS 

maneuver     : 

t  i  me    = 

1?00 

sec 

a 

vdot     = 

0.0 

m/sec /m  i  n 

b 

hdot     = 

0.0 

deq/m  i  n 

c 

fdot     = 

0.0 

hertz/min 

d 

RUOYS   PARAMFFFRS 
First  buoy  :  tvpe  =   DTF4R 


X  -  position  = 
y  -  position  = 


0.0       tm 
-3,0       Wm 


b 
c 


If       OIFAR: 

bearing    error    nnean=  0.0       deq 

St  d    dev    =  S . 0       aeg 


d 
e 


BUOYS   PaRAMFTFRS 

<?nd   buoy  :  tvoe  =   DTFAR 

¥  -  position  =     '4.0 
y  -  position  =    -5,0 


\c  m 


b 

c 


If   DTFAR: 

bearing  error  mp^n- 
st  d  oev  = 


0.0 

s.o 


oeq 
d»q 


d 
e 


Fiaure  E-10  :  List i no  of  problem  variables  ITI 
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SEQUENCE  OF  EVENTS  FRU'^  601  TO  900  SECONDS 


700  sec  -  bean' no  measuremenf  Hy  buoy  nu'nber  1 

700  sec  ••  frequencv  measurernent  by  buoy  nynber  1 

700  sec  -  bearina  measurempor  by  buov  number  2 

700  sec  -  frequencv  neasurempnf-  by  buoy  nui^ber  2 

700  sec  -  time  delay  measurement  bv  buoys  1  and  ? 

700  sec  ""  '''^onte  Carlo  roint  number  7 

^00  sec  -  bearino  measurpment  by  buov  number  1 

000  sec  -  frequency  '"easurement  by  ouoy  number  1 

POO  sec  -  bear i no  measurement  by  buov  number  Z 

POO  sec  -  frequency  measurement  by  buov  number  2 

800  sec  "  time  delay  mpasure'^ent  by  buoys  1  and  2 

POO  sec  -  ^onte  Carlo  coint  number  7 

900  sec  •"  bearinq  measurement  by  buov  number  1 

900  sec  "  frequencv  measurempnt  by  buov  number  1 

900  sec  ~  bear  1  no  r-easurement  by  buov  nu'^bpr  2 

900  sec  -  freauency  ^easurpment  by  buov  number  2 

900  sec  -  time  delay  mpasurement  by  buoys  1  and  ? 

900  sec  ""  f'^onte  Carlo  coint  number  7 


Fiaure  E-11  :  Table  of  events. 
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PWnGPAN'  CHARACTFRISTICS 

Press  the  nun^ber  corresnonHinq  to  the  option  or 

action  desired  and  c/r. 

To  continue  or  return  oress  0  c/r. 


1  •  update  tarcet  everv  second 


?  -  update  tarcet  only  before  'neasurements 


3  -  start  orintinq  on  line  printer 


^  -  stoo  print  inn 


5  -  out  results  of  next  run  on  outout  files 


6  -  modify  parameters  of  the  problem 


7  -  terminate  thp  oroblei 


F inure  E-1?  t  Interruption  table. 
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School/  to  any  interested  reader. 
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